{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" }, "colab": { "name": "503_Stability.ipynb", "provenance": [], "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "Slm9La2lQd_y" }, "source": [ "# Stability of a Multistep method\n", "\n", "#### John S Butler \n", "john.s.butler@tudublin.ie \n", "\n", "[Course Notes](https://johnsbutler.netlify.com/files/Teaching/Numerical_Analysis_for_Differential_Equations.pdf) [Github](https://github.com/john-s-butler-dit/Numerical-Analysis-Python)" ] }, { "cell_type": "markdown", "metadata": { "id": "rjz9aQfcQd_z" }, "source": [ "## Overview\n", "A one-step or multistep method is used to approximate the solution of an initial value problem of the form\n", "\\begin{equation} \\frac{dy}{dt}=f(t,y),\\end{equation}\n", "with the initial condition\n", "\\begin{equation} y(a)=\\alpha.\\end{equation}\n", "The method should only be used if it satisfies the three criteria:\n", "1. that difference equation is consistent with the differential equation;\n", "2. that the numerical solution converges to the exact answer of the differential equation;\n", "3. that the numerical solution is stable.\n", "\n", "In the notebooks in this folder we will illustate examples of consistent and inconsistent, convergent and non-convergent, and stable and unstable methods. \n", "\n", "This notebook focuses on stable and unstable methods. The video below outlines the notebook." ] }, { "cell_type": "code", "metadata": { "id": "tSDjdyrIQd_1", "colab": { "base_uri": "https://localhost:8080/", "height": 336 }, "outputId": "6bdde7a4-1e2b-4398-a499-f3aadf47856c" }, "source": [ "from IPython.display import HTML\n", "HTML('')" ], "execution_count": 1, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "execution_count": 1 } ] }, { "cell_type": "markdown", "metadata": { "id": "BQw_0RUEQd_5" }, "source": [ "## Introduction to unstable\n", "This notebook illustrates an unstable multistep method for numerically approximating an initial value problem\n", "\\begin{equation} \\frac{dy}{dt}=f(t,y), \\end{equation}\n", "with the initial condition\n", "\\begin{equation} y(a)=\\alpha,\\end{equation}\n", "using the Modified Abysmal Kramer-Butler method. The method is named after the great [Cosmo Kramer]( https://en.wikipedia.org/wiki/Cosmo_Kramer) and myself [John Butler](https://johnsbutler.netlify.com).\n", "\n", "## 2-step Abysmal Kramer-Butler Method\n", "\n", "The 2-step Abysmal Kramer-Butler difference equation is given by\n", "\\begin{equation}w_{i+1} = w_{i-1} + h(4f(t_i,w_i)-2f(t_{i-1},w_{i-1})) \\end{equation}\n", "by changing $F$, the Modified Abysmal Butler Method (see convergent and consistent notebooks), the Abysmal Kramer-Butler method is consistent with the differential equation and convergent with the exact solution (see below for proof).\n", "But the most important thing is that method is weakly stable, it fluctuates widely around the exact answer, just like it's name sake Kramer (for examples see any Seinfeld episode)." ] }, { "cell_type": "markdown", "metadata": { "id": "E0xiYTTnQd_6" }, "source": [ "## Definition of Stability\n", "The stability of a numerical method is not as tangable as consistency and convergence but when you see an unstable solution it is obvious.\n", "\n", "\n", "To determine the stabilty of a multistep method we need three definitions:\n", "\n", "\n", "### Definition: Characteristic Equation\n", "Associated with the difference equation \n", "\\begin{equation} w_0=\\alpha \\ \\ \\ w_1=\\alpha_1 \\ \\ \\ ... \\ \\ \\ w_{m-1}=\\alpha_{m-1} \\end{equation}\n", "\\begin{equation}w_{i+1} = a_{m-1}w_{i}+a_{m-2}w_{i-1}+...+a_{0}w_{i+1-m} +hF(t_i,h,w_{i+1},...,w_{i+1-m}),\\end{equation}\n", "is the __characteristic equation__ given by\n", "\\begin{equation}\\lambda^{m} - a_{m-1}\\lambda^{m-1}-a_{m-2}\\lambda^{m-2}-...-a_{0} =0. \\end{equation}\n", "\n", "### Definition: Root Condition \n", "\n", "Let $\\lambda_1,...,\\lambda_m$ denote the roots of the that characteristic equation\n", "\\begin{equation}\\lambda^{m} - a_{m-1}\\lambda^{m-1}-a_{m-2}\\lambda^{m-2}-...-a_{0} =0 \\end{equation}\n", "associated with the multi-step difference method\n", "\\begin{equation} w_0=\\alpha \\ \\ \\ w_1=\\alpha_1 \\ \\ \\ ... \\ \\ \\ w_{m-1}=\\alpha_{m-1} \\end{equation}\n", "\\begin{equation} w_{i+1} = a_{m-1}w_{i}+a_{m-2}w_{i-1}+...+a_{0}w_{i+1-m} +hF(t_i,h,w_{i+1},...,w_{i+1-m}),\\end{equation}\n", "If $|\\lambda_{i}|\\leq 1$ for each $i=1,...,m$ and all roots with absolute value 1\n", "are simple roots then the difference equation is said to satisfy the __root condition__.\n", "\n", "### Definition: Stability\n", "1. Methods that satisfy the root condition and have $\\lambda=1$ as the only root \n", "of the characteristic equation of magnitude one and all other roots are 0 are called __strongly stable__;\n", "2. Methods that satisfy the root condition and have more than one distinct root\n", "with magnitude one are called __weakly stable__;\n", "3. Methods that do not satisfy the root condition are called __unstable__.\n", "\n", "All one step methods, Adams-Bashforth and Adams-Moulton methods are all stongly stable.\n", "\n", "\n", "\n", "## Intial Value Problem\n", "To illustrate stability we will apply the method to a linear intial value problem given by\n", "\\begin{equation}y^{'}=t-y, \\ \\ (0 \\leq t \\leq 2), \\end{equation}\n", "with the initial condition\n", "\\begin{equation}y(0)=1,\\end{equation}\n", "with the exact solution\n", "\\begin{equation}y(t)= 2e^{-t}+t-1.\\end{equation}" ] }, { "cell_type": "markdown", "metadata": { "id": "9UpHqoOlQd_7" }, "source": [ "## Python Libraries" ] }, { "cell_type": "code", "metadata": { "id": "OrZM8puGQd_8" }, "source": [ "import numpy as np\n", "\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "import pandas as pd\n", "\n", "warnings.filterwarnings(\"ignore\")" ], "execution_count": 2, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "ULRF_gMZQd__" }, "source": [ "### Defining the function\n", "$$ f(t,y)=t-y.\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "JendtjxoQd__" }, "source": [ "def myfun_ty(t,y):\n", " return t-y" ], "execution_count": 3, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "C3FvO2V8QeAB" }, "source": [ "## Discrete Interval\n", "Defining the step size $h$ from the interval range $a \\leq t \\leq b$ and number of steps $N$ \n", "\\begin{equation}h=\\frac{b-a}{N}.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=t_0+ih,\\end{equation}\n", "where $t_0=a.$" ] }, { "cell_type": "code", "metadata": { "id": "MG1dErxRQeAC", "colab": { "base_uri": "https://localhost:8080/", "height": 281 }, "outputId": "8683ff93-fe31-4523-e980-c224eb19979a" }, "source": [ "# Start and end of interval\n", "b=2\n", "a=0\n", "# Step size\n", "N=16\n", "h=(b-a)/(N)\n", "t=np.arange(a,b+h,h)\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(t,0*t,'o:',color='red',label='Coarse Mesh')\n", "plt.xlim((0,2))\n", "plt.ylim((-0.1,.1))\n", "\n", "plt.legend()\n", "plt.title('Illustration of discrete time points')\n", "plt.show()" ], "execution_count": 4, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAEICAYAAAA0vXKRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5xVdb3/8ddbQAivgIgoKJiooajJRHZ+KeYVM8XKCn+oWCqppfmrX+doHo8d1LxUR4+/PMfIS6ikGGXS6XgIb1lp5oyh4i0uXhhCHK6KJoJ+fn+s78Bi2Hsu7IFh7Xk/H4/9mPX9ru/6rs/67jV7PrNuWxGBmZmZmRXTVh0dgJmZmZltPCdzZmZmZgXmZM7MzMyswJzMmZmZmRWYkzkzMzOzAnMyZ2ZmZlZgTubMtmCSzpD0h1w5JO3VkTGVI+kmSZd2wHrPlbRI0kpJfVrR/hVJR6Xp70i6edNHuelIGivptx0dR3Mk7Z7eny4dHYtZNXIyZ7aFyScbm6j/n0q6osI+1ksyASLinIi4vLLo2hxHN+DfgGMiYtuIWNKW5SPiexFx1qaJbkOSBqWEvGt7LR8RkyPimPaLsv1FxGvp/Xm/pbaVjpFZZ+RkzszWU7A/ov2AHsBzHR0IFG7szKxKOJkzKyhJj0g6K1dee7RMmeskvSHpTUnPStpf0nhgLPCP6bTXr1P7VyT9k6RngLcldZV0kaS5kt6S9Lykz6a2HwFuAj6R+lie6tc74ifpbElzJC2VNE3Srrl5IekcSbMlLZd0oySV2c7ukq6X9Lf0uj7V7Q28lJotl/RQmeVPk/SqpCWSLmky77uS7kzTPSTdmdotl/SkpH5pXm9Jt6X1L5P0q1R/uKT6NHavA7dJ2io3dksk3SOpd1rlo7l4V0r6ROrnK5JeSH1Pl7RHmbd9g+XLnIo/L43tW5Iul/RhSY+lfeEeSVvn2n9G0sy0zY9JOqDMuhv7vkDSPEmLJX1f0lZp3laS/jmN9RuSbpe0Q5q33tG2tO9eLumPKcbfStqpmW3cS9LvJK1I651SLkazzsjJnFl1OgY4DNgb2AH4IrAkIiYCk4Fr02mvE3LLnAIcD+wYEWuAucChafl/Be6U1D8iXgDOAR5PfezYdOWSjgCuSuvtD7wK3N2k2WeAjwEHpHbHltmWS4BDgIOAA4ERwD9HxF+B/VKbHSPiiBJxDAX+EzgN2BXoAwwos55xaVsHpnbnAH9P8+4Aeqb17Qxcl1tuF6A3sAcwHjgfOAkYmda5DLgxtT0sF++2EfG4pNHAd4DPAX2B3wN3lYlxg+XLtDsWGE42bv8ITAROTdu2P9l7jaSPArcCX03b/GNgmqTuZfoF+CxQAxwMjAa+kurPSK9PAXsC2wI/aqaf/w18mWw8twb+bzPbeDnwW6AX2fv3/5rp16zTcTJnVp1WA9sB+wKKiBciYmELy9wQEfMj4u8AEfHziPhbRHwQEVOA2WSJVGuMBW6NiKciYhVwMdmRvEG5NldHxPKIeA14mCxZK9fXhIh4IyIayBLL01oZx8nAf0XEoymOS4EPyrRdTZbQ7BUR70dEXUS8Kak/cBxwTkQsi4jVEfG73HIfAJdFxKo0ducAl0REfVrnd4GTVf4U7DnAVek9WgN8DziomaNzrXFtRLwZEc8Bs4DfRsS8iFgB3A98NLUbD/w4Ip5I2zwJWEWWBJZzTUQsTe/b9aTEkOx9+re0npVk7/mYZrb7toj4axqzeyj//kP23uwB7BoR70bEH5ppa9bpOJkzq0IR8RDZUZEbgTckTZS0fQuLzc8XJJ2eO/22nOyIzk6lF93ArmRH4xrjWQksAXbLtXk9N/0O2ZGcFvtK07uWaVtq2bXbFRFvpzhKuQOYDtydTqdeq+wGi4HA0ohYVma5hoh4N1feA7g3N24vAO+TXd9Xyh7Av+faLwXE+mPVVoty038vUW4c6z2AbzWuO61/IM2Pb34/yb8Xpd6nrpTf7ta+/5AdXRTwZ0nPSfpKM23NOh0nc2bF9TbZqb9Gu+RnRsQNETEcGEp2uvXbjbPK9Le2Ph0V+gnwdaBPOpU6i+wPanN9NPobWaLQ2N82ZEe9FrSwXIt9AbunutZYSJacNMbRM8WxgXTE7V8jYijwD2SngU8nS156S9rgdHLjok3K84HjImLH3KtHRCwo0bax/VebtP9QRDzWinVVaj5wZZN194yIcqd5ITeerP9elHqf1rB+ItkaG2xjRLweEWdHxK5kp4T/Q1voI3rMOoKTObPimgl8TlLP9IftzMYZkj4m6ePpyNLbwLusO724iOyapuZsQ/ZHtSH192WyI3ONFgED8hfSN3EX8GVJB6Xrr74HPBERr7RlA3N9/bOkvuki+X8B7mzlslOBz0j6ZIp1AmU+9yR9StIwZc9Ce5Ps1N4H6fT0/WQJRC9J3SQdVqqP5CbgysbTpCnu0WleA9n7sGeT9hdL2i+130HSF8r0XWr5SvwEOCftK5K0jaTjJW3XzDLfTuMwEPgG0Hgzwl3A/5E0WNK2ZO/5lHTquC022EZJX5DUeK3jMrJ9s9zpcrNOx8mcWXFdB7xHllhNIruxodH2ZH+ol5Gd7loCfD/NuwUYmk6r/apUxxHxPPBD4PHU/zDgj7kmD5E9DuR1SYtLLP8A2fVpvyA7OvZhYMxGbSVcAdQCzwDPAk+luhala8a+BvwsxbEMqC/TfBey5O9NslOjvyM79QrZNXqrgReBN4ALm1ntvwPTgN9Kegv4E/DxFM87wJXAH9P4HxIR9wLXkJ3efZPsCOhxZbZng+VbHIRmREQtcDbZKfllwByymxiacx9QR/bPxG/I9ifIbqS4g+xu1JfJ/oE4fyNiKrWNHwOekLSSbGy/ERHz2tq3WbVSRHsftTczs2okKYAhETGno2Mxs3V8ZM7MzMyswNolmZM0StJLyh4QelGJ+YdJekrSGkknN5k3TtnDLWdLGperH67sQadzJN0glX6gqJmZmVlnVnEyly4WvpHsGo+hwCnpQZ15r5Fdh/GzJsv2Bi4ju55kBHCZpF5p9n+SXcsxJL1GVRqrmZltvIiQT7GabXna48jcCGBOelDke2RPeR+dbxARr0TEM2x499GxwIz0AMplwAxgVHpI5/YR8afILuq7neyJ6mZmZmaW0x5fCr0b6z9Esp5059ZGLrtbetWXqN+Asu+aHA+wzTbbDN93331buWozMzOzjlNXV7c4IvpW2k97JHMdKn3X5ESAmpqaqK2t7eCIzMzMzFom6dWWW7WsPU6zLmD9J4IPoPVPeS+37ALW/zLstvRpZmZm1mm0RzL3JDAkPfV7a7IHg05r5bLTgWPS08R7AccA09MT19+UdEi6i/V0sgdVmpmZmVlOxclc+qqWr5MlZi8A90TEc5ImSDoR1n61UD3wBeDHkp5Lyy4FLidLCJ8EJqQ6gPOAm8meSD6X7Ot0zMzMzCynqr4BwtfMmZmZZVavXk19fT3vvvtuR4fS6fXo0YMBAwbQrVu39eol1UVETaX9F/4GCDMzM9tQfX092223HYMGDcLP3e84EcGSJUuor69n8ODBm2Qd/jovMzOzKvTuu+/Sp08fJ3IdTBJ9+vTZpEdIncyZmZlVKSdyW4ZN/T44mTMzMzMrMCdzZmZmtkm8/vrrjBkzhg9/+MMMHz6cT3/60/z1r3/t0Jh++tOfIokHHnhgbd2vfvUrJDF16tQ29/fd736XH/zgB+0ZYps5mTMzMzOYPBkGDYKttsp+Tp5cUXcRwWc/+1kOP/xw5s6dS11dHVdddRWLFi2quN8PPmj6Ve9tM2zYMO6+++615bvuuosDDzywoj47kpM5MzOzzm7yZBg/Hl59FSKyn+PHV5TQPfzww3Tr1o1zzjlnbd2BBx7IoYceSkTw7W9/m/33359hw4YxZcoUAFauXMmRRx7JwQcfzLBhw7jvvuz7Al555RX22WcfTj/9dPbff3/mz5/PGWecsXb56667DoC5c+cyatQohg8fzqGHHsqLL75YMrZDDz2UP//5z6xevZqVK1cyZ84cDjrooLXz6+rqGDlyJMOHD+fYY49l4cKFANxwww0MHTqUAw44gDFjxqxt//zzz3P44Yez5557csMNN2z0mG0sP5rEzMysMzj8cDjjjOy1ejUcfTScdRaceipcfDG888767d95By68EMaOhcWL4eST4VvfghNOgNdfh112aXZ1s2bNYvjw4SXn/fKXv2TmzJk8/fTTLF68mI997GMcdthh9O3bl3vvvZftt9+exYsXc8ghh3DiiScCMHv2bCZNmsQhhxxCXV0dCxYsYNasWQAsX74cgPHjx3PTTTcxZMgQnnjiCc477zweeuihDdYviaOOOorp06ezYsUKTjzxRF5++WUgez7f+eefz3333Uffvn2ZMmUKl1xyCbfeeitXX301L7/8Mt27d1+7ToAXX3yRhx9+mLfeeot99tmHc889d4Nnym1KTubMzMw6u/r60vVLlmyS1f3hD3/glFNOoUuXLvTr14+RI0fy5JNPctxxx/Gd73yHRx99lK222ooFCxasPS27xx57cMghhwCw5557Mm/ePM4//3yOP/54jjnmGFauXMljjz3GF77whbXrWbVqVdkYxowZww033MCKFSv44Q9/yPe+9z0AXnrpJWbNmsXRRx8NwPvvv0///v0BOOCAAxg7diwnnXQSJ5100tq+jj/+eLp370737t3ZeeedWbRoEQMGDNhwpZuIkzkzM7PO4JFH1k1367Z+effds1OrTe2+e/Zzp53Wb9/CUTmA/fbbr803FEyePJmGhgbq6uro1q0bgwYNWvt8tm222WZtu169evH0008zffp0brrpJu655x6uv/56dtxxR2bOnNmqdY0YMYJnn32Wnj17svfee6+tjwj2228/Hn/88Q2W+c1vfsOjjz7Kr3/9a6688kqeffZZALp37762TZcuXVizZk2btrtSvmbOzMyss7vySujZc/26nj2z+o10xBFHsGrVKiZOnLi27plnnuH3v/89hx56KFOmTOH999+noaGBRx99lBEjRrBixQp23nlnunXrxsMPP8yrpRJMYPHixXzwwQd8/vOf54orruCpp55i++23Z/Dgwfz85z8HsqTs6aefbjbGq6++eu0RuUb77LMPDQ0Na5O51atX89xzz/HBBx8wf/58PvWpT3HNNdewYsUKVq5cudHj0558ZM7MzKyzGzs2+3nJJfDaa9kRuSuvXFe/ESRx7733cuGFF3LNNdfQo0cPBg0axPXXX88nP/lJHn/8cQ488EAkce2117LLLrswduxYTjjhBIYNG0ZNTQ377rtvyb4XLFjAl7/85bV3tV511VVAdmTv3HPP5YorrmD16tWMGTOm2btUjzvuuA3qtt56a6ZOncoFF1zAihUrWLNmDRdeeCF77703p556KitWrCAiuOCCC9hxxx03enzakyKio2NoNzU1NVFbW9vRYZiZmXW4F154gY985CMdHYYlpd4PSXURUVNp3z7NamZmZlZgTubMzMzMCszJnJmZWZWqpkupimxTvw9O5szMzKpQjx49WLJkiRO6DhYRLFmyhB49emyydfhuVjMzsyo0YMAA6uvraWho6OhQOr0ePXps0ocIt0syJ2kU8O9AF+DmiLi6yfzuwO3AcGAJ8KWIeEXSWODbuaYHAAdHxExJjwD9gb+necdExBvtEa+ZmVm169atG4MHD+7oMGwzqPg0q6QuwI3AccBQ4BRJQ5s0OxNYFhF7AdcB1wBExOSIOCgiDgJOA16OiPyjm8c2znciZ2ZmZrah9rhmbgQwJyLmRcR7wN3A6CZtRgOT0vRU4EhJatLmlLSsmZmZmbVSeyRzuwHzc+X6VFeyTUSsAVYAfZq0+RJwV5O62yTNlHRpieTPzMzMrNPbIu5mlfRx4J2ImJWrHhsRw4BD0+u0MsuOl1QrqdYXeZqZmVln0x7J3AJgYK48INWVbCOpK7AD2Y0QjcbQ5KhcRCxIP98CfkZ2OncDETExImoioqZv374VbIaZmZlZ8bRHMvckMETSYElbkyVm05q0mQaMS9MnAw9FevCNpK2AL5K7Xk5SV0k7peluwGeAWZiZmZnZeip+NElErJH0dWA62aNJbo2I5yRNAGojYhpwC3CHpDnAUrKEr9FhwPyImJer6w5MT4lcF+AB4CeVxmpmZmZWbVRNT4auqamJ2trajg7DzMzMrEWS6iKiptJ+togbIMzMzMxs4ziZMzMzMyswJ3NmZmZmBeZkzszMzKzAnMyZmZmZFZiTOTMzM7MCczJnZmZmVmBO5szMzMwKzMmcmZmZWYE5mTMzMzMrMCdzZmZmZgXmZM7MzMyswJzMmZmZmRWYkzkzMzOzAnMyZ2ZmZlZgTubMzMzMCszJnJmZmVmBOZkzMzMzKzAnc2ZmZmYF1i7JnKRRkl6SNEfSRSXmd5c0Jc1/QtKgVD9I0t8lzUyvm3LLDJf0bFrmBklqj1jNzMzMqknFyZykLsCNwHHAUOAUSUObNDsTWBYRewHXAdfk5s2NiIPS65xc/X8CZwND0mtUpbGamZmZVZv2ODI3ApgTEfMi4j3gbmB0kzajgUlpeipwZHNH2iT1B7aPiD9FRAC3Aye1Q6xmZmZmVaU9krndgPm5cn2qK9kmItYAK4A+ad5gSX+R9DtJh+ba17fQJwCSxkuqlVTb0NBQ2ZaYmZmZFUxH3wCxENg9Ij4KfBP4maTt29JBREyMiJqIqOnbt+8mCdLMzMxsS9UeydwCYGCuPCDVlWwjqSuwA7AkIlZFxBKAiKgD5gJ7p/YDWujTzMzMrNNrj2TuSWCIpMGStgbGANOatJkGjEvTJwMPRURI6ptuoEDSnmQ3OsyLiIXAm5IOSdfWnQ7c1w6xmpmZmVWVrpV2EBFrJH0dmA50AW6NiOckTQBqI2IacAtwh6Q5wFKyhA/gMGCCpNXAB8A5EbE0zTsP+CnwIeD+9DIzMzOzHGU3i1aHmpqaqK2t7egwzMzMzFokqS4iairtp6NvgDAzMzOzCjiZMzMzMyswJ3NmZmZmBeZkzszMzKzAnMyZmZmZFZiTOTMzM7MCczJnZmZmVmBO5szMzMwKzMmcmZmZWYE5mTMzMzMrMCdzZmZmZgXmZM7MzMyswJzMmZmZmRWYkzkzMzOzAnMyZ2ZmZlZgTubMzMzMCszJnJmZmVmBOZkzMzMzK7B2SeYkjZL0kqQ5ki4qMb+7pClp/hOSBqX6oyXVSXo2/Twit8wjqc+Z6bVze8RqZmZmVk26VtqBpC7AjcDRQD3wpKRpEfF8rtmZwLKI2EvSGOAa4EvAYuCEiPibpP2B6cBuueXGRkRtpTGamZmZVav2ODI3ApgTEfMi4j3gbmB0kzajgUlpeipwpCRFxF8i4m+p/jngQ5K6t0NMZmZmZp1CeyRzuwHzc+V61j+6tl6biFgDrAD6NGnzeeCpiFiVq7stnWK9VJJKrVzSeEm1kmobGhoq2Q4zMzOzwtkiboCQtB/Zqdev5qrHRsQw4ND0Oq3UshExMSJqIqKmb9++mz5YMzMzsy1IeyRzC4CBufKAVFeyjaSuwA7AklQeANwLnB4RcxsXiIgF6edbwM/ITueamZmZWU57JHNPAkMkDZa0NTAGmNakzTRgXJo+GXgoIkLSjsBvgIsi4o+NjSV1lbRTmu4GfAaY1Q6xmpmZmVWVipO5dA3c18nuRH0BuCcinpM0QdKJqdktQB9Jc4BvAo2PL/k6sBfwL00eQdIdmC7pGWAm2ZG9n1Qaq5mZmVm1UUR0dAztpqamJmpr/SQTMzMz2/JJqouImkr72SJugDAzMzOzjeNkzszMzKzAnMyZmZmZFZiTOTMzM7MCczJnZmZmVmBO5szMzMwKzMmcmZmZWYE5mTMzMzMrMCdzZmZmZgXmZM7MzMyswJzMmZmZmRWYkzkzMzOzAnMyZ2ZmZlZgTubMzMzMCszJnJmZmVmBOZkzMzMzKzAnc2ZmZmYF5mTOzMzMrMDaJZmTNErSS5LmSLqoxPzukqak+U9IGpSbd3Gqf0nSsa3ts6S6Ohg0CCZPrnyjNtbkyVkMW23lWBxLMeJwLI7FsVRXLFtKHI6lxViGw/B26S8iKnoBXYC5wJ7A1sDTwNAmbc4DbkrTY4ApaXpoat8dGJz66dKaPku9hkMERPTsGXHnnbHZ3Xlntu7GOByLY9nS43AsjsWxVFcsW0ocjqVVsQyHiArzsMi2puJk7hPA9Fz5YuDiJm2mA59I012BxYCatm1s15o+S72G59+k3XePGDky4o47ssF7++2sfPfdWXn58qz8i19k5YaGrDxtWlZeuDAr339/Vn7ttaw8Y0ZWnjs3Kz/ySFZ+8cWI7t3X31EaX3vsEfGXv2Tt//KXrP2f/5yVn302K//xj1n5xRez8iOPZOW5c7PyjBlZ+bXXsvL992flhQuz8rRpWbmhISvvtFPpWAYOzObfcUfW/r33svJtt2XlRhMnRhx55LryjTdGjBq1rnz99REnnLCu/P3vR3zuc+vKV10V8aUvZdN77FE6lm22Wdf+oosizj57Xflb34o477x15W98I3s1Ou+8rE2js8/O+mh0xhkRl166rjx2bMSECeVj+dCHsm1odMIJ2TY2GjUqG4NGRx6ZjVGjkSOzMYzIxrSlfa+5fWVj9r2RI7N9KCLbp0aOzPaxiJb3vX79Ssey667Z/Lbue7/4RVZevjwr3313Vn777azc3L5X7v3p0WPdWLdl34vI3vexY9eVL7002z8aldv3ysWy7bYbt+81+tKXshgbfe5zLe97vXqV31/auu9V+rnX3L7b1n2v0s+9vn2b/5xry74XUdnnXrn9pWfPde1bu+812pjPvR12KP/+bMy+V8nnXnP7Snv8zW3L594uuzT/OVfp39y2fO7l9pX2Suba4zTrbsD8XLk+1ZVsExFrgBVAn2aWbU2fAEgaL6lWUu16M+bPL9V801q1qnT9a69t3jgAliwpXV9fv3njgPLb//bbmzcOKB/L3/++eePYkvaVN94oXb9w4eaNA8pv/7vvbt44oHwsK1du3jgAli0rXd8R+8uWtO8uXly6fkv6nHvnnc0bx4oVpes7+76yaFHp+i3pc64SlWaDwMnAzbnyacCPmrSZBQzIlecCOwE/Ak7N1d+S+muxz1Kv4U0z/82t3H9mjsWxbKlxOBbH4liqK5YtJQ7H0qpYtqQjcwuAgbnygFRXso2krsAOwJJmlm1Nn+X17AlXXtnq5u3myiuzdTsWx1KUOByLY3Es1RXLlhKHY2lbLJWqNBskuwZuHtkNDI03K+zXpM3XWP8GiHvS9H6sfwPEPLKbH1rss9RreGOW3REXNDa6884sBsmxOJZixOFYHItjqa5YtpQ4HEuLsbTXkTlFllRVRNKngetTInZrRFwpaQJQGxHTJPUA7gA+CiwFxkTEvLTsJcBXgDXAhRFxf7k+W4qjpqYmamtrW2pmZmZm1uEk1UVETcX9tEcyt6VwMmdmZmZF0V7JnL8BwszMzKzAnMyZmZmZFZiTOTMzM7MCczJnZmZmVmBO5szMzMwKzMmcmZmZWYE5mTMzMzMrMCdzZmZmZgXmZM7MzMyswJzMmZmZmRWYkzkzMzOzAnMyZ2ZmZlZgTubMzMzMCszJnJmZmVmBOZkzMzMzKzAnc2ZmZmYF5mTOzMzMrMCczJmZmZkVWEXJnKTekmZImp1+9irTblxqM1vSuFTXU9JvJL0o6TlJV+fanyGpQdLM9DqrkjjNzMzMqlWlR+YuAh6MiCHAg6m8Hkm9gcuAjwMjgMtySd8PImJf4KPA/5J0XG7RKRFxUHrdXGGcZmZmZlWp0mRuNDApTU8CTirR5lhgRkQsjYhlwAxgVES8ExEPA0TEe8BTwIAK4zEzMzPrVCpN5vpFxMI0/TrQr0Sb3YD5uXJ9qltL0o7ACWRH9xp9XtIzkqZKGlguAEnjJdVKqm1oaNiojTAzMzMrqhaTOUkPSJpV4jU63y4iAoi2BiCpK3AXcENEzEvVvwYGRcQBZEfyJpVbPiImRkRNRNT07du3ras3MzMzK7SuLTWIiKPKzZO0SFL/iFgoqT/wRolmC4DDc+UBwCO58kRgdkRcn1vnktz8m4FrW4rTzMzMrDOq9DTrNGBcmh4H3FeizXTgGEm90o0Px6Q6JF0B7ABcmF8gJYaNTgReqDBOMzMzs6pUaTJ3NXC0pNnAUamMpBpJNwNExFLgcuDJ9JoQEUslDQAuAYYCTzV5BMkF6XElTwMXAGdUGKeZmZlZVVJ2qVt1qKmpidra2o4Ow8zMzKxFkuoioqbSfvwNEGZmZmYF5mTOzMzMrMCczJmZmZkVmJM5MzMzswJzMmdmZmZWYE7mzMzMzArMyZyZmZlZgTmZMzMzMyswJ3NmZmZmBeZkzszMzKzAnMyZmZmZFZiTOTMzM7MCczJnZmZmVmBO5szMzMwKzMmcmZmZWYE5mTMzMzMrMCdzZmZmZgXmZM7MzMyswCpK5iT1ljRD0uz0s1eZduNSm9mSxuXqH5H0kqSZ6bVzqu8uaYqkOZKekDSokjjNzMzMqlWlR+YuAh6MiCHAg6m8Hkm9gcuAjwMjgMuaJH1jI+Kg9Hoj1Z0JLIuIvYDrgGsqjNPMzMysKlWazI0GJqXpScBJJdocC8yIiKURsQyYAYxqQ79TgSMlqcJYzczMzKpOpclcv4hYmKZfB/qVaLMbMD9Xrk91jW5Lp1gvzSVsa5eJiDXACqBPqQAkjZdUK6m2oaGhgk0xMzMzK56uLTWQ9ACwS4lZl+QLERGSoo3rHxsRCyRtB/wCOA24vS0dRMREYCJATU1NW9dvZmZmVmgtJnMRcVS5eZIWSeofEQsl9QfeKNFsAXB4rjwAeCT1vSD9fEvSz8iuqbs9LTMQqJfUFdgBWNKaDTIzMzPrTCo9zToNaLw7dRxwX4k204FjJPVKNz4cA0yX1FXSTgCSugGfAWaV6Pdk4KGI8FE3MzMzsyZaPDLXgquBeySdCbwKfBFAUg1wTkScFRFLJV0OPJmWmZDqtiFL6roBXYAHgJ+kNrcAd0iaAywFxlQYp5mZmVlVUjUd8KqpqYna2tqODsPMzMysRZLqIqKm0n78DRBmZmZmBavYCXQAAAorSURBVOZkzszMzKzAnMyZmZmZFZiTOTMzM7MCczJnZmZmVmBO5szMzMwKzMmcmZmZWYE5mTMzMzMrMCdzZmZmZgXmZM7MzMyswJzMmZmZmRWYkzkzMzOzAnMyZ2ZmZlZgTubMzMzMCszJnJmZmVmBOZkzMzMzKzAnc2ZmZmYF5mTOzMzMrMAqSuYk9ZY0Q9Ls9LNXmXbjUpvZksaluu0kzcy9Fku6Ps07Q1JDbt5ZlcRpZmZmVq0qPTJ3EfBgRAwBHkzl9UjqDVwGfBwYAVwmqVdEvBURBzW+gFeBX+YWnZKbf3OFcZqZmZlVpUqTudHApDQ9CTipRJtjgRkRsTQilgEzgFH5BpL2BnYGfl9hPGZmZmadSqXJXL+IWJimXwf6lWizGzA/V65PdXljyI7ERa7u85KekTRV0sAK4zQzMzOrSl1baiDpAWCXErMuyRciIiRFiXatMQY4LVf+NXBXRKyS9FWyo35HlIlvPDAeYPfdd9/I1ZuZmZkVU4vJXEQcVW6epEWS+kfEQkn9gTdKNFsAHJ4rDwAeyfVxINA1Iupy61ySa38zcG0z8U0EJgLU1NRsbDJpZmZmVkiVnmadBoxL0+OA+0q0mQ4cI6lXutv1mFTX6BTgrvwCKTFsdCLwQoVxmpmZmVWlFo/MteBq4B5JZ5LdjfpFAEk1wDkRcVZELJV0OfBkWmZCRCzN9fFF4NNN+r1A0onAGmApcEaFcZqZmZlVJa1/z0Gx1dTURG1tbUeHYWZmZtYiSXURUVNpP/4GCDMzM7MCczJnZmZmVmBO5szMzMwKzMmcmZmZWYE5mTMzMzMrMCdzZmZmZgXmZM7MzMyswJzMmZmZmRWYkzkzMzOzAnMyZ2ZmZlZgTubMzMzMCszJnJmZmVmBOZkzMzMzKzAnc2ZmZmYF5mTOzMzMrMCczJmZmZkVmJM5MzMzswJzMmdmZmZWYE7mzMzMzAqsomROUm9JMyTNTj97lWn3P5KWS/qvJvWDJT0haY6kKZK2TvXdU3lOmj+okjjNzMzMqlWlR+YuAh6MiCHAg6lcyveB00rUXwNcFxF7AcuAM1P9mcCyVH9damdmZmZmTVSazI0GJqXpScBJpRpFxIPAW/k6SQKOAKaWWD7f71TgyNTezMzMzHK6Vrh8v4hYmKZfB/q1Ydk+wPKIWJPK9cBuaXo3YD5ARKyRtCK1X9y0E0njgfGpuErSrLZtQqewEyXGzjwuJXhMSvO4lOZxKc3jsiGPSWn7tEcnLSZzkh4Adikx65J8ISJCUrRHUG0REROBiQCSaiOiZnPHsKXzuJTmcdmQx6Q0j0tpHpfSPC4b8piUJqm2PfppMZmLiKOaCWKRpP4RsVBSf+CNNqx7CbCjpK7p6NwAYEGatwAYCNRL6grskNqbmZmZWU6l18xNA8al6XHAfa1dMCICeBg4ucTy+X5PBh5K7c3MzMwsp9Jk7mrgaEmzgaNSGUk1km5ubCTp98DPyW5kqJd0bJr1T8A3Jc0huybullR/C9An1X+T8nfJNjWxwu2pVh6X0jwuG/KYlOZxKc3jUprHZUMek9LaZVzkA15mZmZmxeVvgDAzMzMrMCdzZmZmZgVWmGRO0ihJL6Wv+NrgGrrmvgJM0sWp/qXc9XqF14ox+aak5yU9I+lBSXvk5r0vaWZ6Tdu8kW9arRiXMyQ15Lb/rNy8cenr6WZLGtd02SJrxbhclxuTv0panptXlfuLpFslvVHu+ZTK3JDG7BlJB+fmVfO+0tK4jE3j8aykxyQdmJv3Sqqf2V6PXdhStGJcDpe0Ive78i+5ec3+/hVVK8bk27nxmJU+S3qnedW8rwyU9HD6G/ycpG+UaNN+ny8RscW/gC7AXGBPYGvgaWBokzbnATel6THAlDQ9NLXvDgxO/XTp6G3aTGPyKaBnmj63cUxSeWVHb0MHjssZwI9KLNsbmJd+9krTvTp6mzbXuDRpfz5wayfYXw4DDgZmlZn/aeB+QMAhwBPVvq+0clz+oXF7geMaxyWVXwF26uht6KBxORz4rxL1bfr9K9KrpTFp0vYEsqdTdIZ9pT9wcJreDvhrib9F7fb5UpQjcyOAORExLyLeA+4m+8qvvHJfATYauDsiVkXEy8Cc1F/RtTgmEfFwRLyTin8ie5ZftWvNvlLOscCMiFgaEcuAGcCoTRTn5tbWcTkFuGuzRNaBIuJRYGkzTUYDt0fmT2TPxuxPde8rLY5LRDyWths6z2dLa/aXcir5XNqitXFMOsXnCkBELIyIp9L0W8ALrPuWq0bt9vlSlGRu7dd7Jfmv/tqgTWQPIW78CrDWLFtEbd2uM8n+A2jUQ1KtpD9JKvmdugXV2nH5fDqsPVXSwDYuW0St3rZ0On4w8FCuulr3l5aUG7dq3lfaqulnSwC/lVSn7OsWO5tPSHpa0v2S9kt1nX5/kdSTLCH5Ra66U+wryi77+ijwRJNZ7fb5Uul3s1oBSDoVqAFG5qr3iIgFkvYEHpL0bETM7ZgIN7tfA3dFxCpJXyU7ontEB8e0JRkDTI2I93N1nXl/sTIkfYosmftkrvqTaV/ZGZgh6cV09KYzeIrsd2WlpE8DvwKGdHBMW4oTgD9GRP4oXtXvK5K2JUtgL4yINzfVeopyZK7x670a5b/6a4M2Wv8rwFqzbBG1arskHUX2PbonRsSqxvqIWJB+zgMeIfuvoRq0OC4RsSQ3FjcDw1u7bIG1ZdvG0ORUSBXvLy0pN27VvK+0iqQDyH5/RkfE2q9bzO0rbwD3Uh2XtbRKRLwZESvT9H8D3STthPcXaP5zpSr3FUndyBK5yRHxyxJN2u/zpaMvEmzlhYRdyS4AHMy6i0f3a9Lma6x/A8Q9aXo/1r8BYh7VcQNEa8bko2QX3Q5pUt8L6J6mdwJmUz0X47ZmXPrnpj8L/ClN9wZeTuPTK0337uht2lzjktrtS3ZRsjrD/pK2aRDlL2g/nvUvUP5zte8rrRyX3cmuP/6HJvXbANvlph8DRnX0tmzGcdml8XeHLDF5Le07rfr9K+qruTFJ83cgu65um86yr6T3/Xbg+mbatNvnSyFOs0bEGklfB6aT3RV0a0Q8J2kCUBsR08i+AuwOZV8BtpQsoSO1uwd4HlgDfC3WP31USK0ck+8D2wI/z+4F4bWIOBH4CPBjSR+QHZ29OiKe75ANaWetHJcLJJ1Itj8sJbu7lYhYKuly4MnU3YRY/5RAYbVyXCD7vbk70idKUrX7i6S7yO5A3ElSPXAZ0A0gIm4C/pvsjrM5wDvAl9O8qt1XoFXj8i9k1yT/R/psWRMRNUA/4N5U1xX4WUT8z2bfgE2kFeNyMnCupDXA34Ex6Xep5O9fB2xCu2vFmED2T/NvI+Lt3KJVva8A/ws4DXhW0sxU9x2yf4Ta/fPFX+dlZmZmVmBFuWbOzMzMzEpwMmdmZmZWYE7mzMzMzArMyZyZmZlZgTmZMzMzMyswJ3NmZmZmBeZkzszMzKzA/j9WPAs8R7U6kgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "1V_764XPQeAE" }, "source": [ "## 2-step Abysmal Kramer-Butler Method\n", "\n", "For this initial value problem 2-step Abysmal Kramer-Butler difference equation is\n", "\\begin{equation}w_{i+1} = w_{i-1} + h(4(t_i- w_i)-2(t_{i-1}-w_{i-1})) \\end{equation}\n", "by changing $F$, the Modified Abysmal Butler Method, is consistent and convergent.\n", "\n", "For $i=0$ the system of difference equation is:\n", "\\begin{equation}w_{1} = w_{-1} + h(4(t_0-w_0)-2(t_{-1}-w_{-1})) \\end{equation}\n", "this is not solvable as $w_{-1}$ is unknown.\n", "\n", "For $i=1$ the difference equation is:\n", "\\begin{equation}w_{2} = w_{0} + h(4(t_1-w_1)-2(t_{0}-w_{0})) \\end{equation}\n", "this is not solvable as $w_{1}$ is unknown. $w_1$ can be approximated using a one step method. Here, as the exact solution is known,\n", "\\begin{equation}w_1=2e^{-t_1}+t_1-1.\\end{equation}\n" ] }, { "cell_type": "code", "metadata": { "id": "duQPu0rrQeAE" }, "source": [ "### Initial conditions\n", "IC=1\n", "w=np.zeros(len(t))\n", "y=(2)*np.exp(-t)+t-1\n", "w[0]=IC\n", "w[1]=y[1]" ], "execution_count": 5, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "mS-2OWTOQeAG" }, "source": [ "### Loop" ] }, { "cell_type": "code", "metadata": { "id": "s5h0LwfOQeAG" }, "source": [ "for i in range (1,N):\n", " w[i+1]=(w[i-1]+h*(4*myfun_ty(t[i],w[i])-2*myfun_ty(t[i-1],w[i-1]))) \n" ], "execution_count": 6, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "OGxetHJ1QeAI" }, "source": [ "### Plotting solution" ] }, { "cell_type": "code", "metadata": { "id": "33atL94wQeAI" }, "source": [ "def plotting(t,w,y):\n", " \n", " fig = plt.figure(figsize=(10,4))\n", " plt.plot(t,w,'^:',color='red',label='Abysmal Kramer-Butler (N)')\n", " plt.plot(t,y, 'o-',color='black',label='Exact')\n", " plt.xlabel('time')\n", " plt.legend()\n", " plt.title(' Abysmal Kramer-Butler')\n", " plt.show " ], "execution_count": 7, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "KrDezp2bQeAK" }, "source": [ "The plot below shows the Abysmal Kramer-Butler approximation $w_i$ (red) and the exact solution $y(t_i)$ (black) of the intial value problem. \n", "\n", "The numerically solution initially approximates the exact solution $(t<.5)$ reasonably but then the instability creeps in such that the numerical approximation starts to widely oscilate around the exact solution." ] }, { "cell_type": "code", "metadata": { "id": "s46QiSk4QeAK", "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "outputId": "100d6c8f-ebc7-4736-817a-a69b5c09d610" }, "source": [ "plotting(t,w,y)" ], "execution_count": 8, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAEWCAYAAABysAOLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXhTZdoG8PvpTgtlK1S2toCyihQpqCgiAi7gvoGCAo7yKa7jCLgwbojjqKjj4IzigiwVEWZAxgEVxxUBZUdZZWkLpWwFCrTQtM3z/fEmJW3TNm2TnKS9f9eVqzknJ+c8b9I0T99VVBVEREREVFKI1QEQERERBSImSURERERuMEkiIiIicoNJEhEREZEbTJKIiIiI3GCSREREROQGkySiICci34nIPVbHUR4RGSUiy6yOI1iISJqIDLQ6DiJikkQUFMTYJSKbrY7Fm0QkSURURMIc2yIifxeRrSLSyur4KiMiz4lIgYicdNy2iMjNVXj+RyLyoi9jJKLqY5JEFBwuBdAcQDsR6WV1ML4gIiEA3gVwGYB+qprp5pgwf8flcm1xxFjaXFWtr6r1ATwKYLaIxPs5PEtfG6LaikkSUXAYCeAzAIsd90trLyK/iMhxEflMRJoAgIj8V0Qecj1QRDaKyI2OL/03ROSg43m/isi5jmM+EpF/iMgSRw3JTyJyloi8KSJHHTU9PVzO+YSI7BSREyKyWURurGL5QgFMB5AC4DJVPeA47yjHtd8QkWwAz4lIexH5RkSyReSwiKSKSCOXWNJEZJyjnLki8oGIxDvKckJEvhaRxi7HXygiy0XkmIhsEJHLXB77TkQmi8hPAPIAtKuoEKr6JYATANq7xF+iqdFRc3a2iIwBMBzAeMdr/J/S5xOREJfXNltEPnV5b521cH8QkQwA31TlBSeiyjFJIgpwIhIN4BYAqY7bMBGJKHXYXQDuBtACQCGAtxz7ZwAY4XKu7gBaAfgvgCtgaqg6AGgI4DYA2S7nvA3ARABxAPIBrACw1rE9H8DrLsfuBNDXcZ7nYWpTWlShmKkAOgK4XFWzSz12AYBdAOIBTAYgAP4CoCWAzgDaAHiu1HNuBjDIUbZrASwB8BSAZjB/9x52vB7O1+JFAE0APA7gXyLSzOVcdwIYA6ABgPTyCuBIOocAiABQabOoqk5zlPsVR03UtW4OewjADQD6Ocp7FMDbpY7pB/M6XFnZNYmoapgkEQW+m2CSlK9gvtDDAQwpdcwsVf1NVXMB/BnAbSISCmARgA4ico7juDthmodsAApgvvg7ARBV3aKqWS7nXKCqa1T1NIAFAE6r6kxVLQIwF0BxTZKqzlPVfapqV9W5AH4H0LsKZbwCwDxVPebmsX2q+ndVLVTVU6q6Q1WXqmq+qh6CSdb6lXrO31X1gKPJ7kcAP6vqOpeyOGMfAWCxqi52xL4UwGoAg13O9ZGqbnJcv8BNfLeJyDEAJ2Fe75fKKUd13AfgaVXdq6r5MMngLaWa1p5T1VxVPeWlaxKRA5MkosA3EsCnji/p0wD+hbJNbntc7qfDJFJxjuPnAhjh6E9zO4BZAKCq3wCYClMzcVBEpolIrMt5DrjcP+Vmu75zQ0TuEpH1jiarYwDOhalx8tQ1AJ4VkbvdPOZaNjiazj4RkUwROQ5gtptreRp7IoBbnXE7Yr8EpkauzPVdOmifFJEEx+5PVbWRqsbANLPdJSL/51GpK5cIYIFLbFsAFMHUqpWJj4i8i0kSUQATkdYALodJcvaLyH6YprfBIuKaGLRxuZ8AU0t02LE9A6bvywAAeaq6wnmgqr6lqj0BdIFpmhpXjRgTAbwH4EEATVW1EYDfYJrFPLUcplnsbyJyR6nHtNT2S4593VQ1FqY2qCrXcrUHphaukcstRlVfdnd9Zwdtxy2j9MlUNQ2mac/ZdJYLINr5uIicVUnZ3MV3dan4okp1aq/sHERUTUySiALbnQC2w/TXSXbcOgDYC1Mr5DRCRLo4+i+9AGC+o1kMjqTIDmAKHLVIACAivUTkAhEJh/kyP+04rqpiYL6oDznOOxqmJqlKVPV7mKbFaVLxMPoGME1bOY4+RVVO7FzMBnCtiFwpIqEiEiUilzmS0ypzPO8qAJscuzYA6CoiySIShbJ9pw6g4s7g7wCY7EhEISLNROT66sRGRFXHJIkosI0E8A9V3e96g/nydG1ymwXgIwD7AUTB0THZxUwA3WCSAqdYmBqgozBNdNkAXq1qgKq6GSYBWwHzpd8NwE9VPY/jXEsBDAUwQ0TcdWQGTMfw8wHkwPTR+nd1ruW43h4A18N06j4EU3MzDlX72zjU2QQHYBVM2Z93nH87TNL6NUw/rdKTan4AoIujOW2hm3P/Daaf01cicgLASpiO7ETkB6LKmlqi2k5E7gIwRlUvsToWIqJgwZokolrO0QQ3FsA0q2MhIgomTJKIajERuRKmGekAgI8tDoeIKKiwuY2IiIjIDdYkEREREbnhkwUR4+LiNCkpyRenJiIiIvKqNWvWHFbVZqX3+yRJSkpKwurVq31xaiIiIiKvEhG36zKyuY2IiIjIDSZJRERERG4wSSIiIiJywyd9ktwpKCjA3r17cfr0aX9dkshjUVFRaN26NcLDw60OhYiIAoTfkqS9e/eiQYMGSEpKgkh1F+wm8j5VRXZ2Nvbu3Yu2bdtaHQ4REQUIvzW3nT59Gk2bNmWCRAFHRNC0aVPWchIRBZKsLKBfP2D/fstC8GufJCZIFKj4u0lEFGAmTQKWLTM/LcKO20RERBRYsrKA6dMBu938tKg2qc4lSQsXLoSIYOvWrcX7vvvuO1xzzTUWRmXUr1+/0v2LFy9Ghw4dkJ7udt4rv7vsssvQsWNHJCcno3Pnzpg2rfKF5hcuXIjNmzcXb48aNQrz58+vURwLFy7ECy+8AAB47rnnEB0djYMHDxY/7nwNbTYbLr30UhQWFtboekRE5EOTJpkECQCKiiyrTQrsJMkH7ZFz5szBJZdcgjlz5njtnP7yv//9Dw8//DCWLFmCxMTEEo/560vf3XVSU1Oxfv16/PTTT5gwYQJsNluF5yidJHkjhldeeQVjx44t3o6Li8OUKVPKHBcREYEBAwZg7ty51b4+ERH5kLMWyfldYrNZVpsU2EmSl9sjT548iWXLluGDDz7AJ598UuKx48ePY8iQIejYsSPuu+8+2O12fPjhh3j00UeLj3nvvffwxz/+Ebm5uRgyZAi6d++Oc889t/gLNykpCU8++SSSk5ORkpKCtWvX4sorr0T79u3xzjvvFMcwYMAAnH/++ejWrRs+++wzj2L/4YcfcO+99+Lzzz9H+/btAZgamPvuuw8XXHABxo8fj19++QUXXXQRevTogT59+mDbtm0AgI8++gg33HADBg0ahKSkJEydOhWvv/46evTogQsvvBBHjhwBAOzcuRNXXXUVevbsib59+xbXtpW+TkWvb0xMDEJDQwGUrAGbP38+Ro0aheXLl2PRokUYN24ckpOTsXPnzhLnWLNmDfr164eePXviyiuvRFZWFgBTY/Xoo48iJSUFf/vb30o8Z/v27YiMjERcXFzxvrvvvhtz584tLpurG264AampqR697kRE5GeutUhOVtUmqarXbz179tTSNm/eXHJHv36q06eb+zab2Z41y2zn5qpeeKFqeLgqoBoVpXrRRar/+pd5/NAhc/yiRWY7K6vM9dyZPXu23n333aqqetFFF+nq1atVVfXbb7/VyMhI3blzpxYWFurAgQN13rx5euLECW3Xrp3abLbi52zcuFHnz5+v99xzT/F5jx07pqqqiYmJ+o9//ENVVR999FHt1q2bHj9+XA8ePKjNmzdXVdWCggLNyclxFOOQtm/fXu12u6qqxsTEuI07LCxMGzdurBs2bCixf+TIkTpkyBAtLCxUVdWcnBwtKChQVdWlS5fqTTfdpKqq06dP1/bt2xfHEhsbq//85z+L43zjjTdUVfXyyy/X7du3q6rqypUrtX///m6v46pfv37aoUMH7datm0ZFRek777xT/JhreebNm6cjR44sPt+8efNKlGPevHlqs9n0oosu0oMHD6qq6ieffKKjR48uvs7999/v9vX58MMP9bHHHivefvbZZ/XVV1/V559/Xp955pkysRQWFmpcXFyZ85T5HSUiIv9LTjbf/aVvyck+uySA1eomn/HbPElVlp5uXhbAZJRe6IMzZ84cPPLIIwCAYcOGYc6cOejZsycAoHfv3mjXrh0A4Pbbb8eyZctwyy234PLLL8fnn3+Ozp07o6CgAN26dUNkZCT+9Kc/YcKECbjmmmvQt2/f4mtcd911AIBu3brh5MmTaNCgARo0aIDIyEgcO3YMMTExeOqpp/DDDz8gJCQEmZmZOHDgAM4666xy4w4PD0efPn3wwQcflKlFufXWW4trbnJycjBy5Ej8/vvvEBEUFBQUH9e/f//iWBo2bIhrr722OM6NGzfi5MmTWL58OW699dbi5+Tn57u9TmmpqalISUnBoUOH0KdPH1x11VVlmgM9sW3bNvz2228YNGgQAKCoqAgtWrQofnzo0KFun5eVlYVmzcos3oyHH34YycnJePzxx0vsDw0NRUREBE6cOIEGDRpUOU4iIvKhdeusjqCYdUnSd9+duR8eXnI7Jwc4ehRw9j2x2cx2nz5mOy6u5PEVJBhOR44cwTfffINff/0VIoKioiKICF599VUAZYeAO7fvuecevPTSS+jUqRNGjx4NAOjQoQPWrl2LxYsXY+LEiRgwYACeeeYZAEBkZCQAICQkpPi+c7uwsBCpqak4dOgQ1qxZg/DwcCQlJVU6P09ISAg+/fRTDBgwAC+99BKeeuqp4sdiYmKK7//5z39G//79sWDBAqSlpeGyyy4rfqx0LK5xFhYWwm63o1GjRli/fr3bGJzXefvtt/Hee+8BMJ3IXTVr1gznn38+fv75ZyQmJpZ4TT2Zg0hV0bVrV6xYsaLCGEqrV68ecnJyyuxv1KgR7rjjDrz99ttlHsvPz0dUVFSlMRERkQXsdiDE+h5B1kfgjg/aI+fPn48777wT6enpSEtLw549e9C2bVv8+OOPAIBffvkFu3fvht1ux9y5c3HJJZcAAC644ALs2bMHH3/8MW6//XYAwL59+xAdHY0RI0Zg3LhxWLt2rcdx5OTkoHnz5ggPD8e3337r8Si16Oho/Pe//0Vqaio++OCDcs/dqlUrAKYfUlXExsaibdu2mDdvHgCTsGzYsKHMcQ888ADWr1+P9evXo2XLliUey8vLw7p164r7TMXHx2PLli2w2+1YsGBB8XENGjTAiRMnypy7Y8eOOHToUHGSVFBQgE2bNlUae+fOnbFjxw63jz322GN49913S3T2zs7ORlxcHJcgISIKVOPGAR06nGlRskhgJkkrVpzp1e5kswHLl1f7lHPmzMGNN95YYt/NN99cPMqtV69eePDBB9G5c2e0bdu2xLG33XYbLr74YjRu3BgA8Ouvv6J3795ITk7G888/j4kTJ3ocx/Dhw7F69Wp069YNM2fORKdOnTx+bpMmTfDFF1/gxRdfxKJFi8o8Pn78eDz55JPo0aNHtUa7OROw7t27o2vXrh53Kh8+fDiSk5PRs2dPjBo1qrgJ8+WXX8Y111yDPn36lGg2GzZsGF599VX06NGjRMftiIgIzJ8/HxMmTED37t2RnJyM5R6855deeinWrVsHdfNhiouLw4033lii6fDbb7/FkCFDPCobERFZ4MILgaFDAYsn+hV3Xyw1lZKSoqtXry6xb8uWLejcubPXr+UP11xzDf74xz9iwIABVodC5XjkkUdw7bXXYuDAgZUee9NNN+Hll19Ghw4dSuwP5t9RIiKqPhFZo6oppfcHZk1SgDh27Bg6dOiAevXqMUEKcE899RTy8vIqPc5ms+GGG24okyAREVGAKCwEPPh77g+BO7otADRq1Ajbt2+3OgzyQHx8fPHIwopERETgrrvu8kNERERULRs3AikpwH/+A1jcNYI1SURERBQ4mjYFnnsO6NbN6khYk0REREQBJDERcEyrYzXWJBEREVHgyMwsO8LdIkySiIiIKHAMHAgMG2Z1FADqWHNbaGgourm0cQ4bNgxPPPGEV869fv167Nu3D4MHD/bK+YiIiOqkZ58FmjSxOgoAAVyTlJqaiqSkJISEhCApKckrq7bXq1eveLbo9evXey1BAkySVHqZDiIiIqqiYcOAK66wOgoAAZokpaamYsyYMUhPT4eqIj09HWPGjPFKolRaTk4OOnbsiG3btgEwi9s61ya7//77kZKSgq5du+LZZ58tfs6qVavQp08fdO/eHb1790ZOTg6eeeYZzJ07F8nJyZg7d67X4yQiIqr19u8Hdu4suzSZRSxpbnv00UfLXUgVAFauXFliGQnArAv2hz/8oTiBKS05ORlvvvlmhdc9deoUkpOTi7effPJJDB06FFOnTsWoUaPwyCOP4OjRo7j33nsBAJMnT0aTJk1QVFSEAQMGYOPGjejUqROGDh2KuXPnolevXjh+/Diio6PxwgsvYPXq1Zg6daqnLwMRERG5+uADYOJE4MQJoH59q6MJzD5JpROkyvZ7ytncVtqgQYMwb948PPDAAyUWdf30008xbdo0FBYWIisrC5s3b4aIoEWLFujVqxcAszAsERERecEttwBt2wZEggRYlCRVVuOTlJSE9PT0MvsTExPx3XffeT0eu92OLVu2IDo6GkePHkXr1q2xe/duvPbaa1i1ahUaN26MUaNG4fTp016/NhERETl07GhuASIg+yRNnjwZ0dHRJfZFR0dj8uTJPrneG2+8gc6dO+Pjjz/G6NGjUVBQgOPHjyMmJgYNGzbEgQMHsGTJEgBAx44dkZWVhVWrVgEATpw4gcLCQjRo0AAnTpzwSXxERES1nt0OLF0KZGdbHUmxgEyShg8fjmnTpiExMREigsTEREybNg3Dhw+v0XmdfZKctyeeeALbtm3D+++/jylTpqBv37649NJL8eKLL6J79+7o0aMHOnXqhDvuuAMXX3wxALP219y5c/HQQw+he/fuGDRoEE6fPo3+/ftj8+bN7LhNRERUHXv2mFFt8+dbHUkxUVWvnzQlJUVXr15dYt+WLVvQuXNnr1+LyFv4O0pEZKFTp4CVK4EOHYBWrfx6aRFZo6oppfcHZMdtIiIiqmPq1QP697c6ihICsrmNiIiI6pjvvzc1SQHErzVJqgoR8ecliTzii2ZnIiKqgokTzc8ff7Q2Dhd+S5KioqKQnZ2Npk2bMlGigKKqyM7ORlRUlNWhEBHVXXPmAMePWx1FCX5Lklq3bo29e/fi0KFD/rokkceioqLQunVrq8MgIqq7AvBvsN+SpPDwcLRt29ZflyMiIqJgsXMn8MMPwI03Ao0aWR1NMXbcJiIiImstXQrcfXfANbdVmiSJSJSI/CIiG0Rkk4g874/AiIiIqI64917g998DrsnNk+a2fACXq+pJEQkHsExElqhqYI3TIyIiouAUGgqcfbbVUZRRaU2SGicdm+GOG8dLExERkXf89a/A8uVWR1GGR32SRCRURNYDOAhgqar+7OaYMSKyWkRWcwQbEREReSQvz8yR9N13VkdShkej21S1CECyiDQCsEBEzlXV30odMw3ANMCs3eb1SImIiKj2iY4GTp4ECgqsjqSMKo1uU9VjAL4FcJVvwiEiIqI6JzISqF/f6ijK8GR0WzNHDRJEpB6AQQC2+jowIiIiqgPmzjV9kgKQJzVJLQB8KyIbAayC6ZP0uW/DIiIiojrh66+BWbOsjsKtSvskqepGAD38EAsRERHVNe+9F5D9kQDOuE1ERERWCw+3OgK3mCQRERGRNXbuBEaPBjZvtjoSt5gkERERkTX27gW+/BI4fdrqSNzyaJ4kIiIiIq/r1w/Yt8/qKMrFmiQiIiIiN5gkERERkTXuuQeYOtXqKMrFJImIiIiskZEBZGdbHUW52CeJiIiIrPHVV1ZHUCHWJBERERG5wSSJiIiI/O/994FBg4D8fKsjKReTJCIiIvK/kBBzi4y0OpJyMUkiIiIi/7v7bjORZABjkkRERETkBpMkIiIi8q/9+4H27YFFi6yOpEJMkoiIiMi/8vOBXr2A5s2tjqRCnCeJiIiI/CsxEfjkE6ujqBRrkoiIiMi/VK2OwCNMkoiIiMi/Bg8Ghg61OopKsbmNiIiI/Ovyy4GoKKujqBSTJCIiIvKvceOsjsAjbG4jIiIKVFlZQL9+Zsh8bWGzAUVFVkfhESZJREREgWrSJGDZMvOztkhNBerXB9LTrY6kUkySiIiIAlFWFjB9OmC3m5+1pTapa1fgwQeBVq2sjqRSTJKIiIgC0cSJJkECTPNUbalN6t0bePVVICzwu0UzSSIiIgo0WVnAzJmm/w5gftaW2qT9+zlPEhEREVXTpEmASMl9taE26fRp08wWJOUI/LouIiKiumbFCqCgoOQ+mw1YvtyaeLzFbgemTjVNbkGASRIREVGgmTEDCAkxnZxL1ygFs+ho4P77rY7CY2xuIyIiCjTPPw9ceaXpu/Pee8A111gdkXdkZACHDlkdhcdYk0RERBRo/v53YOdOU5tUVGSa3k6fDoqlPCo0bhywerUpWxAQ9UEP85SUFF29erXXz0tERERBbPly4OBB4IYbrI6kBBFZo6oppfezJomIiCiQfPQR0Lw5MHhwyf2qwd8/qU8fqyOoEvZJIiIiCiSvvmoSJVfjxwOXXWZFNN5z7BiwciVw6pTVkXiMNUlERESBZP16k1C4Ovts0y8pmGuTfvwRuO460+R20UVWR+MRJklERESBJDwcaNas5L4xY6yJxZsuugj47DPgvPOsjsRjbG4jIiIKFMOHA5984v4xVeDkSf/G401xcaYmKSbG6kg8xiSJiIgoEJw4AWzfbkZ/uXPNNSbJCFZffGHKF0TY3EZERBQIGjQAVq0qf/HXoUPPLHgbbFRN/CNGAG+/bXU0HmOSREREFAicnbLL65h9113+jcfbli8PuskwK21uE5E2IvKtiGwWkU0i8og/AiMiIqozDh8GWrYEFiyo+Ljjx83SHsFGxKxD17691ZFUiSd9kgoB/ElVuwC4EMADItLFt2ERERHVISdOAP37A23bVnxcr17AI0FYV7FqFTB/vlliJYhU2tymqlkAshz3T4jIFgCtAGz2cWxERER1Q9u2wMcfV37cyy8DTZv6Ph5vmz7dlO/mm62OpEqqtHabiCQB+AHAuap6vNRjYwCMAYCEhISe6enp3ouSiIiotsrPB3JyzFIktVVuLrB3L9Cxo9WRuFXe2m0eTwEgIvUB/AvAo6UTJABQ1WmqmqKqKc1KT4JFRERE7i1dCpx1FrBiReXHqgLr1gGbNvk+Lm+KiQnYBKkiHiVJIhIOkyClquq/fRsSERFRHdK1K/D888D553t2/BVXAK+95tuYvOn4ceCVV4AdO6yOpMoq7ZMkIgLgAwBbVPV134dERERUh7RtC/z5z54dK2I6QLdr59uYvGnLFmDCBKBLF7MGXRDxZJ6kiwHcCeBXEVnv2PeUqi72XVhERER1wO7dwP79wAUXACEe9oDp18+3MXnbBRcAR48CkZFWR1JlnoxuWwYgSJccJiIiCmDvvWeaog4dAho39uw5+fnAwoXAOed43kRntUaNrI6gWrh2GxERkVXGjweWLPE8QQJMk9vo0cDs2b6Ly5umTgXmzLE6imrhsiRERERWadQIGDSoas+JiDAj3IJl9uoPPwQ6dABuv93qSKqMSRIREZEVFi0yfXXuuqv89drKE0zD6desMU2EQYjNbURERFaYMQN4442qJ0gAcOSImTZgzRrvx+VtIkG3sK0TkyQiIiIrzJ8PfPll9Z4bFga8+KJnE1Ba6dtvgYcfNjVmQYhJEhERkRVEgPj46j03NtYkHg8+6N2YvG3rVmDWLCA62upIqoVJEhERkb89/rhpaquJ+vW9E4sv3X+/aRoMwjmSACZJREQUyLKyzOSJ+/dbHYn3qALbtwNpaTU7z65dwPDhwIYNXgnLZ6rT5ypAMEkiIqLANWkSsGyZ+VlbiJiRbW++WbPzREWZPj8ZGd6Jy9sKC4HrrgMWB+8CHUySiIgoMGVlAdOnA3a7+VlbapPsdvOzpjUsLVsCmZnAtdfWPCZfyM4G0tOBnByrI6k2JklERBSYXnjhTEJRVFQ7apPsdjOxYk1rkZwCuSkrPt40BQbhJJJOTJKIiCjwZGWZmZptNrNts9WO2qSTJ4GBA703W/aKFUBKiunjRF7HJImIiAKPu1qj2lCbFBsLvPOO95rImjQBGjQwyVegefpp4L77rI6iRrgsCRERBZ4VK87UIjnZbMDy5dbE4y179gBt2njvfB07ms7bgaiw0NyCGGuSiIgosGRkAJ06mYRC9cwtMzN4Vr535/ffgYQEYOZM75/b2XcrkPz1r8D771sdRY0wSSIiosCydi3w9dclv/hVzXxJgT7DdEUaNwamTAH69/fueefPN81uBw5497zEJImIiALMDTeY2qSEhDP7RIB33wU++MC6uGoqLg547DHvNrcBwDnnAMOGlW2etNK33wLJycCWLVZHUiNMkoiIKDAUFgIrV5r79eqVffzyy4F27fwbk7ccPQp89ZVvEpnu3U1ncG8nXzURFmbmcWra1OpIaoRJEhERBYbp04GLLgJ+/rn8Y/btA0aPBjZv9l9c3vDZZ8CVVwK//uq7axw+7LtzV1Xfvmam7ebNrY6kRpgkERFRYBg+3MyN1Lt3+cdERACffw6sW+e/uLxh6FCTNJx/vm/O/8orwFlnBeZUAEGMUwAQEZH1VIHoaFNLVJG4ODPqLSrKP3F5S716wNVX++78AwcC4eGBM8rt7LOBUaOAiROtjqRGWJNERETWWrsWuPBCYMcOz453JkiHDvkuJm9asQL4+9+BvDzfXeP884E//tFMVmm1ggKTEHbubHUkNcYkiYiIrJWdbb5YmzXz/Dmvv246cR854ru4vOU//zE1KmE+brzJywuMZsjwcJMU3nyz1ZHUGJMkIiKy1qBBwJo1QMOGnj/niiuAceN8n3h4w0svmbXVIiJ8e50JE0yH6YIC316nMlZf34uYJBERkTVOnABSU00/mqquZn/uudQBILcAACAASURBVMAzzwRG85In4uN9f4177gE+/dT316nM+PFAYqLpZxbkmCQREZE1PvoIGDGi+sPiVYHvvzej3QLV66+bCST9kTB07w4MHmyau6zUr59J2Kqa+LpITU1FUlISQkJCkJSUhNTUVC8G6LkgqKckIqJa6YEHgB49zJd7dT31lJmE8pprvBeXN2VkALt21ShhqJLt24Hdu82cTFa54QZzq6bU1FSMGTMGeY6O7unp6RgzZgwAYPjw4V4J0VOiPshuU1JSdPXq1V4/LxFRnZSVZZadmDvXzIVTG5w+7Z1h/Lt3m9fE3QzdgULVf0nSXXcBX34J7N/vv2u6KiwEcnOr1r+slMTERGRkZLjdn5aWVoPgyicia1Q1pcx+JklERAFu7Fizbtl99wFvv211NDX33/8C//d/wNKl3hsm7vwusyIxKI/dDoT4uVfL9u3mmu3bW/NarFkDpKSYGcavu87tIaqKQ4cOYdeuXW5ve/bscfs8EYHdR/NAlZcksbmNiCiQZWWZRV3tdrNsx5//HPy1Sc2bA5dcYr7IvSE9HbjlFuCFF3w7YWNVDRoEdO0KvPWW/67ZoYP/ruVOfDzw8ss43bkz0rZuLTcRys3NLfG0li1bol27dujfvz8+++wz5OTklDl1guuCx37CJImIKJA988yZRVGLioBJk4K/NqlXL+CTT7x3vhYtgAYNzOsTKFRNORMT/X/tJUuAnBzTRFsDqampePrpp5GRkYGEhARMnjy5uE+QquLgwYPlJkGZTz4J15aq6OhotGvXDu3atcOAAQOK77dr1w5JSUmo59JcWrpPkvP5kydPrlF5qoPNbUREgchuN/1K2rc3/Xec6tUzHYGDsTZp3z5TK/b444HdhyjYXX+96au1cWO1T5Gamop7770Xp06dKt4XFhaG8847DzabDbt27SqRxABAq1at0K5lS7Tr0MHcXBKh+Ph4SBWa/ypK0HyBfZKIiIKFzQbccYdpRtq48UxNktPNNwPz51sTW0384x9mOPymTd5ranNVWAgsXw5ceqn3z11VGRlAmzbW9Avavx9o0qTSyStPnz6N9PR0pKWlFd92796NtLQ0rFq1ym3/n/DwcFx99dVo164d2rdvX6I2KCoqyvQx69QJWLDAV6XzCfZJIiIKFuHhZnTQvn1lEySg+vMKWW3sWDNU31d9S95808zCvWWL+aK2yqlT5vqPPQa8+KL/r++oZczPz0dGRobbJCgtLQ1ZWVklnhYeHo6EhAQkJSWV20G6sLAQn332WfnXfumlGo1sCzRMkoiIAkVmpllmIz4eeP9997UQrsPJvTWM3teKioADB4CWLX2XIAFm1fn27a3vvKwK/O1vQM+eNTpNZU1ONpsNe/bscZsApf32G/YdP16iX1BoaGhxEnT11VcjKSmp+Na2bVu0aNECoaGhAICkpCSkp6eXianSztM33lijMgcaNrcREQWCwkLgvPNMgvTNN5U308yYYWopfvjBdFwOZNOmmRXqV6+uFSvD+4O7zsvh4eHo3bs3QkJCkJaWhszMzBI1PiEhIWjTpo1JfLZuRdu4OCSNG1ecCLVq1QphHq51V17n6WnTppXfNygzEzh2zNSiOZKtYMHmNiKiQBYWZpawiI/3rB9L585mtupgaNoYNMh01vZXE9i775rRXePH++d6roqKzBxBgwaZEXflUFUcPXoUGRkZxbf09PTi+6tWrUJRqdF6BQUFWLlyJfr06YPLLrsMbdu2LVEb1Lp1a4Q7lySx2Wq0oK4zEapS5+mPPgImTjRr8tWvX+1rBxLWJBERWWnRIjOSrQbLOCA/36y8Xku+mGrsrrtM5+Uvv6zx+mFVHmG1YgXQpw8KZs7E3ksuKZEElU6GSs8VFBERgYSEBCQkJOCbb75xe3pfTqhYY7t2AevWmYEFQabao9tE5EMA1wA4qKrnenIxXyZJqampeHrCBGRkZiKhdWtMfvllv6/lQkTkFXa7GYkVGgp89131vtBVTWfo/Hzgq6/8P8NzRX79FZgyBXj1VaBZM/9d1wt9tSpqbrrjjjtw7NixMrU/GRkZyEhPR8bOndh3+DBKf782a9asOAlyvSUmJiIhIQHNmjVDiOP9K69PkMdLc9jtZpHZXr2A+++v0WtRF9Skue0jAFMBzPR2UFVVZtG7vXstW/SOiPysNq1fpmqaZcLCTNNMZGT1azxEgNtvNzVJgZQgAaYP0ldfmUTJn5wJUm6uSUCrmDDl5uZi/PjxZeYBysvLw6hRo3Dffffh5MmTJR6LiIhAmzZtkJiYiEFDhpRJhNq0aYPo6GiPY5g8eXLNJlQMCTHTEPhrMku73Sw307u3aTKuJTxqbhORJACfW12TVFFve3f7iagWqS3rl6masuTkALNnez+x2b3bfDEGSsJ06pQ1E0dmZpqO8M8/Dzz4IADTDygnJwd79+5FZmYm9u7dW+Lm3Hf06NEKT/3II4+USYKaN2+OkLQ0M5P4PfeYpVdqyN8TKtZIejqQlAS8845Zly/I1GgySU+SJBEZA2AMACQkJPT0RdISEhJSpvrSqVd0NLr06YMugwahyznnoGuTJki8+GKEeNiTn4gCVFYW8OyzwKxZphklmGecdnr5ZeD4cWDyZO9ONrh3L9CtG/DwwyY5sEpenpmryMdD4F3Z7XYcPny4ZPIzaxb2Nm6MzPz84n2l+wEBQHx8PFq3bo3WrVujVatWaN26NaZMmYLs7Owyx1bY3DVjBjB69JlEtS7Jzzf9kRITA3+0pRs+T5Jc+awmqXVrpGdmltnfIDoaF0RFYXNICPYdPly8v15EBDqfey66JCSgS14eugwdii59+6Jdu3bFc0EQUYAqKjJNJb//bkZFhYSYYfKhoaaG4OefzaSLweLIEeDgQd+O8FIFXnsNuO02a7+kX3rJLMS7bRtw9tnVOoW7PkFRUVEYM2YM2rVrV6Y2KDMzE7ZSE2+GhoaiZcuWZRIg1+2WLVsiws0osGoNgQfMe+yFWiSvyM0FBg4E7rzT1F5SuWrFFACTO3XCmMxMuLYSRwP4Z58+GL50KQDg6NGj2LJsGTbPn4/N9ephc1oavl+2DLMPHzZt4wAiw8PRKTQUXa64Al169UKXFi3QpUULtL/8coRX0nbNjuNEPqZq+h41amSa1+rXN0OZneuXFRUBGzYAhw+b/1hdJ1cMZMOGATt3Alu3+i65EzEzTgPmdVm1yvQR8bcHHjCTRlaSINlsNhw4cABZWVnYt28fsrKyim+pqakl1g0DzDIab731FgAgMjKyONnp06eP2wQoPj4eodnZQGoq8OijVfo9qdYQeCBwEiQAiIkBWrXyzzQR//uf+b0OhCVhvElVK70BSALwmyfHqip69uypPpGcrLMBTQRUHD9nA6rJyZU+NWfPHl3500/64Ycf6uO33qqDmzXTpMREBVB8Cw8P13PPPVdv69tXn+vXTz/9+GP97bffND8/X1VVZ8+erdHR0SWeEx0drbNnz/ZNeYnqCrtd9Zdfzmw/9ZTqSy+Z+/ffrxoRoWq+9s0tIkJ17FjVY8dUzztPdckSa+Kuig0bVL/5xn/Xmz3bvFY1vObs2bM1MTFRRUQTExMr/3tnt2teXp7u3LlTly1bpvPmzdO33npLn3zySR05cqReccUV2q1bN42Liyvxt9R5ExGNj493+5jz8UOHDqndbvesADNmqIqorlpVo9ehUp9+qnrXXarHj/v2OoHq4otV+/a1OopqA7Ba3eQznkwBMAfAZQDiABwA8KyqflDRc4JpnqTc3Fxs/fprbPrqK2xu0ACbN2/G5mXLsOvoUThfmbCQEJwTEYE0kTL/2QBAYqtWSFu2zHRaA8x/vGFh5uZttWmED1VNbX7v33oLeOQRYPPmsjMy9+gBrF9f9jnJycC8eWYpirfeAs4/33z2ajJSzNt+/NH003j4Yf9f22YDPvjAdKKtZifu8pq8/vSnP6Fr164lan6ysrKQtWMH9mVmIqfUJIiAWUH+rLPOQosWLYpvLVu2LLHdokULNG/eHGFhYTUfAu9UUGD6sHXsWJ2XwHN/+xswfbp5vwPl98/Jbjc1sL5snj540AxGOOcc313Dh2rUJ6mqgilJKs+p48exbdcubN68GZvmzMHmrVuxcMeOco+/smFDJA0bhsTERCTNno3EsDAkLVmCs846CyE33ABER5tRD4AZ+dC0KfDXv5rtp54yzQYPPWS2//53U0V6001m+1//MmseXXSRaVd+5x3zxfDhh757ASjw1JbRXYBJZt5/3zQF9e5tms4WLzZD2Wvyh3zcOJOY/PBDjWYb9ppRo0zfqbVrrRnh5ZSdbRLQvn1L7Lbb7Thy5AgOHjxY4nbgwAEcPHgQs2bNcvuPoauoqKgzSY7djha7d6PF2LFo0aZNiSSoadOmxXMAeaLafYIq4uum2UBs+t2zx/yj8dpr5veR3CovSfKo+ayqN581t1kssVUrt9W/9SIiNKVDB7fVxxEREXp206Y64Jxz9A9/+INOmjRJZ15+uf4wYoSmp6drYWGh6hVXqN5335kLnX226ogRZ7Zbt1YdPVpnT51asqnx7bfN42PGqP7rX2eOP3DANF/URvv2qV56qWpWltWR+J7dfuZ9TE01TQaAar16qps2qR45Ym181eEsT26ualyc6rhx3j3/7NmqTzxxZjsz07vn99Tp02d+Hj5co1NVtbnr5MmTumvXLl25cqUuWrRI33//fZ3crZs+Ehmpt99yiw4YMEC7deum8fHxGhIS4vZvWkhIiDZv3rzCJq9NmzbpkSNHyjZ7efFvT5Wb+iryxBOqt93mtdhKKCryzXm9oajINFn/+KPvrrF9u+o//6mane27a/gYqtvcVh21oSbJndSBAzHmf/8r03F82sCBxR3Hc3NzkZ6ejrS0NLc/9+/fX+KcYWFhaN26NZKSkkwtVFISEhMSkJSQgMR27dCmTRuEZ2Uh9fPPMeaRR5BXWHjm2mFhmPb++xg+ebKpnRo/3vyHHh0NvPCCWUPHZjPDp2+91TRHON/vQPtvx1MjR5q5ZWpDbUppp06Z9yc62ixtcP31wH/+A1xwgalVXLjQPB4RAXTvbpqgDh8GYmPNyKnYWN808XrLa6+ZwRPOpSIyM02Nqa/s2mWa7qZOBe6913fXKW3SJGDJEuDrr817WQPualMiIiJw6623ok2bNmVqgA4ePFhmAkSnBtHRaO5ozoqPj0fz5s3LvTVp0gShoaGeN3kdPGhq8G66KXD/trz8spnL5+23vT+H1L33mt/nxYu9e95gMW2aadZNTzcd9oNQrRjdZrXhjjkzngaQASABwGQAw12mHYiJiUGXLl3QpUsXt+c4ffo0MjIySiRPzvtff/019u3bV2IuqJCQELRq1QoHDxxAvkuCBAB5hYUY/8QTuP733xETEwMBTNvzW28BF15oDsrMNF9OnTubJGn3bjN8euZM8wctO9t8+V59tWnSq4w3+sW4Vklv3w4cOHCmGWDWLDNL7NNPm+377zfb//2vufbs2aaM06ebIcbffGNGQQ0eXL1YrLR3rxnO3qKFGfXUqZPpQ3LXXUDbtqZM9eqZci9ZcibBtdmAjRtNIhwba/Y99JBp0tmyxWzn5pqRLVbbvh1o396UMzbWjPw5dcokD75MkAAgLg6YMAEYMsRs79ljksgqzuFS5Qn9unY11yrV3FdYWIijR48iOzsbR44cKfGzvPt79+4tMzeczWZDamoqwsLCSiQ2HTp0KDf5adasGeo5m/u++AJo1w7o0KHSsns86/PUqSYJ2b79TN/MQPPEE747d7dugT830KFD5m9CDRN3t+6913yHtG7t/XNbjDVJAcZms2HPnj1laqFmzqx4VZh69eqV+x9ifFwcmsfFoXmrVmheUIC4GTMQdu+95oP9v/+ZeTS++Qbo3x9YvtzU0syaZWor9u83Q5Z790bqggV4+r77kHHyJBLq18fkd97B8GHDgGPHTB8rwKzVtGEDMGKE2U5NNTUHzvj/9CeTYO3da7bvucckAM75r+65x3zZr11rtt980/yX+tJLpk/Oe++ZuXIiIsyx331nOmT++9/m+JtuAlJSTD8vwNQmtG4dGP1TVq82yUKPHiZRiI01tX+TJ5vE79lnzcKQycklnzd2rEmeXOeAcZbfWZv2+efmvbrnHrOdkmIS41mzzHZhof9rmX76ySS/8+YFxoKXw4YB335rku7ISI+eUl6/mDfffBMDBgw4k9RkZSH7t9+Q3aBBuYlPTk5OudcJDQ1FkyZN0KRJEzRt2hRNmzZFkyZNMGPGDLfHiwiKioogVa21OXXKDMvv3RtYsMDj16DSJLGwEPjlF6BPn6rFY4Vt28yQ+No2+KEiq1ebNdwWLKjZQsq1GDtuB7mkiAikFxSU2d80NBTjX3qp3I6XhaVqnwDzB7Zp06YmiWrWDM3r10d8QgKat2yJ5rm5aP7114h/+mk079YNzb/6CvXHjsXHU6ZgzMSJyHPpxBldrx6mDRmC4fPnnxlV9PzzwHPPmRElYWGmFmv27DMjPhYuNPedswFv3QqcOGE+wBXJyjL//TrnygFMLcvWreanc/HM2283s/w+/ripeWnY8MzoJ1XglVdMUljDmYA9snSpqdFx/lFq187U5s2fb7ZTU8126dFcpVU0umvdurL7Vc1Im7POMolBYSHQpo3p1PzYY2eO8XaziKops81mFlwtKgJef900kdZw7hivLM+wcyewYQNODx6MnJwc5Hz6KY736IGc/Hyz7bgdP368+P6cOXPKbb4qT6NGjcokO6Xvl94XGxvrtlOz10Z4udq0yfzj4I25c+x281n3MOm0XHa2qTF/4AHzu+kNe/aYcwbyBMX5+eYfzltuMTW73vbCC0C/fuYWpJgkBbnqjPRQVRw7dqxE0uQukXLejh075vY8URERKCgqQpGbYb2NGjTAS1ddhdirr0ZsXBxi7XbEhoYitmNHxDZqhNjYWER64w/o2LFInTYNTxcVnWnqDA3F8P/7v/L7JhUWmlqrs882/XoOHzZf1q+/biaWy8kxNR0vvghcd535cs/IME1d7v7gVdbU+NlnprnLWa1/1VWmKdGZyKxcaZIVXzczlZaTY5LSq68GBg0ytXYXXGBq5a6+2qNTeJSkqJpm3qgo4PvvvRZ+ecPQJ0yYgAsvvNBtclPevuPHj5eZldmdmJgYNGzYEPv27Sv3mI8++uhMwmO3o2l6OhoNHYowL9bY+WSEl1NhoamhfeSR6idMM2aYPljffRc8TS3z5pkvc29N+tipE3DuuWf++alr8vJMS8IzzwBPPml1NNXGJKkW8PVih/n5+Th06JDbROq1116r9nkjIiIQGxtbo9uS3r3x0N69ZTvNJyRgeFXWCczLMzUcDRqYhOj++02tU//+pomvZ0/g009NR/c9e0xtz4gRSP3++7JNjdHR5tiPPza1Mg8/bDpa79pltvfsMX88vNAHwKvv/e+/mz5fzzxj/rh//71p9ps50+1cMhV+UTdtCjz/PPTrr1EQHo68rVuRW78+8oqKkJeXh9zcXOTl5VX5vuu+jIwM2O12j4omIoiNjUXDhg2Lb6W3i/ft24eGycloGB+PhmlpiM3NRcPbbkNsw4bFiU6FNTmffWbe+5df9mlnZZ997letAi6+2CTLI0dW7xzffWf6B06fHjgL6vqT3W7+BjRpAlxxhdXRVCw/3zS7nX++96ejKCoy/2RaOc1FDTFJohpJatAA6SdPltnfpn59/PL77zhx4gSOHz9e7dtp12a0KggPD0f37t0RERGByMhIREZGlrhfervC+6dPI3LDBkRceCEiW7ZE5M8/I3LCBHz18MN47t13cTo/v/i69aKiMOX663HLxo3AokVAw4aQ06dNLYpImb4irttVvT9v3jw89NBDJearqVevHv7yl79gyJAhKCgogM1mK/HTk/vF+7Ztg+2bb1Bw880oEIHt119RsHMnCi6+GDZVLFy40G2TU1hYGOJiY5GXk4NcwG1NY2Wio6MRExOD6Ojocu/PcvarKkVEsGzZshLJT0xMTJXm4ik2YoSZW2n7dvMeOlSYIKalmTnLVq0K3v4tO3ZUe221oPbbb2aeunfeCYwBDv6weLEZxPDtt8Bll1kdTcBhkkQ1kpqYiDEZGTWvySmHzWarMNG67777yn3u4MGDkZ+fD5vNhvz8/DL3S29TWaGhoQgPD0dERATCi4oQbrMhokULhIeHY+fOneU+79577jEJTf36FSY67u5HRUV51PHYJ/1ySissBNLSTMJgt5tawXvvBbp3L1uT8+yzGD56tGlezM42o+iC3fbtwD/+AUyZ4lnfmq1bzWCBhx8OjEERVbVsmZliY/Fi0/RcXf/5j2lidvaJDGTHj58ZoOPNtdxSU02y/eyz3junBZgkUY35urmvIt76olRVFBYWlptAuUuubrzhBpT3KZk6dWrxeV2vUfqaNbn/+OOPl1uemTNnmsQmPPxMkuPy09P7FdW+lFeL6NUkpQI+7ZfjzrZtZpTW3/8O3HGH2efsj3bddSaZ+Omn4K09cmfqVNP59uefTZ+8yrzwgkmofv89sBZ0rYq8vJo1hR86BMTHm2br557zWlhBZ+xYMyra3eCSIMIZtymoWbm4cGL9+m5nHU6sX9/n11ZVTSy1EHPx9RMT/XL92bNna3RYWMnXPizMrws7e3XmZU8cO3ZmFuUZM1S7dlUNCVG95RbV4cNV8/J8e30rHDpUtePT030Th79Vd5Zou1117VrVPXu8G48vZWaqfvihakGBd88byDOOewjlzLjNJImCht+/KJ3XTUjQ6FIJSjSgsxMS/HN9CxNEVVXdt09nh4WVXBInPLxuLA2jqvr44yZBci4JU9vL/dZbqn/5i/vHTp1SPXjQv/H40t13q557bu1dxqm01FTze7x2rdWRBJzykqQ6OByBgtXw4cORlpYGu92OtLQ0vzX1DU9Px7TZs5GYmAgRQWJiIqbNnu2VvlgeXX/4cEybNq3k9X3V1OTOpEkYHhKCNAB2AGkAhouYod91QW7umYk4i4pqd7lVTZPbzz+bvlmlTZliZuquYGqEoHLddWY5jaoOOsjLMyNEK+ivF5AGDzYLHXfv7p3z7dhh5qb77TfvnC8AsU8SEVWsqpNZ1iblTWK6a1ft6pPkqqDATGkQFlZ20tEtW8yszc4Z7euqZcvMCLEvvjCT09ZVy5aZZZT+/e+yKwUEmfL6JLEmiYgqtm6d+bIsfavtCRJgao1K16jU9tqk8HCTIB0/biYf/fxzkyz26wc0blz7EiS73SR+Vfl9vuQS03E7GGeYXrPGTC7rjQqSSy4x/zAEeYJUESZJRETlWbGi5Jp5gNlevtyaePwtLw84edJMuvrDD8DEiVZH5H2nTpk1D//5z6o9r3Fjk1AGm19+MetF1pYmUx9jcxsREblnt5uldRISzFxStbWpcdMmM9u8J0vK/PyzmYjy9deBpCSfh+Z1ublmdnRvzI59551mBu8//rHm57IYm9uIiKhqQkJM06JzHq3a2tTYteuZPliV2bfP9NFr3Nj3cflCTIz3lg/JyTG1jbUYa5KIiMi9utRxfcUK0+y2ZImpOatI6Q7twWbhQmDp0vIXB6+DWJNERERVU5c6rrdqBcTGmqVmyuN8LYI5QQLMsjJffWX6Y1GFmCQREZF7danjekKCKW+PHuUf89prwHnnmX49wWz8eLOkTE2a3d59F+jbt2QtYy3kQS81IiKqk+rCNA+lnTplJlzs2bPsY23bmgVtY2L8H5c3VbBWo8fq1QMaNQKiomp+rgDGPklEREROd95p+iVlZNRsAdxAN2UKsHYtkJpqdSQBgX2SiIiIKjNuHDB/ftmmqH37gPx8a2LyBZvN1Jr5oKKkNmGSRERE5HTeeWbJkdKds++/330TXLB68kmznEh1OqHv3w+0bGmeX8uxTxIREZGr/HzglVfM/Ek33WT2PfhgxSPfglV1pjMoKACuuAJo3do3MQUQJklERESuwsOBefOAI0fOJEmDBlkbky/84Q/A0aNVrxFq0wb46COfhBRomCQRERG5Cgkx0wE4R7F9/TVw9tnBuQxJRTp1MmvzVVVRERAa6v14AhCTJCIiotKcCdLmzcDgwea2cKG1MXnbuHHVe95VVwH16wMLFng3ngDEJImIiMid778HBgwwM23X1ukAVM3kmPXre/6ca68FIiJ8F1MAYZJERETkjrNjsqqpRdq/v/atWde7N9C+PfDJJ54/5+GHfRdPgOEUAERERO5MmXKm701tXbNu9Gjg+us9P/706bJL1dRiTJKIiIhKy8oCpk8/kxDYbGZ7/35r4/K2sWOB22/3/Pg5c0x/rbQ0n4UUSJgkERERlTZpkumL5Kq21iYdPWpmFPdE9+7AhAl1Yo4kgH2SiIiIylqxomyzks0GLF9uTTy+YrebPkm33gq8+27lx59/vrnVEUySiIiISlu3zuoI/CMkBHj7bTMPlCfS000tUh2ZJ4nNbURERHXZ7bcDvXpVftzp00C7drWzybEcrEkiIiKqy4qKgDVrgCZNKq5RUgXeew/o0cN/sVmMNUlERER1WUEB0LcvMG1axcfVqwfcfTeTpNJE5CoR2SYiO0TkCV8HRURERH4SFQX897/AY49VfNyuXcCePf6JKUBUmiSJSCiAtwFcDaALgNtFpIuvAyMiIiI/GTiw8tnEn3wSuOwyv4QTKDzpk9QbwA5V3QUAIvIJgOsBbPZlYEREROQnJ0+apVd69QI6dnR/zLhxwMGD/o3LYp40t7UC4Fq/ttexrwQRGSMiq0Vk9aFDh7wVHxEREfmazQbceSewYEH5x6SkAIMH+y+mAOC10W2qOg3ANABISUlRb52XiIiIfKxJE2DLFuCcc9w/fvQosGGDSZTq1/dvbBbypCYpE0Abl+3Wjn1ERERUW3TqVP4kkcuWAf37A7/+6t+YLOZJkrQKCsrfYQAAChtJREFUwDki0lZEIgAMA7DIt2ERERGRX+3bB0ycCGzfXvaxSy4BvvwSOO88/8dloUqTJFUtBPAggC8BbAHwqapu8nVgRERE5Ec2G/Dyy8DatWUfa9wYuOIKICbG/3FZyKM+Saq6GMBiH8dCREREVklMBI4dc9/naNEisxBu167+j8tCnHGbiIiIABH3CZIqMGIE8M9/+j8mizFJIiIiImPjRuC224CMjJL7164FHn/cmpgsxAVuiYiIyFAFfv7ZJEkJCWafSMUL39ZiTJKIiIjIOO88ID295L6VK4GdO4Fhw8qfIqCWYnMbERERGSJl982cCTz4IBBS91KGuldiIiIiKt8XXwDJycCRI2b7zTdNnyR3CVQtxySJiIiIzoiNBeLjgexssx0RAbRta21MFmGSRERERGf06WNm1z7nHCAnB3jhBWDbNqujsgSTJCIiIirLbgd+/x147jlgxw6ro7EEkyQiIiIq6d13gWbNzAzbJ08CAwdaHZElmCQRERFRSZ07A3fcAeTmAtHRQGSk1RFZgvMkERERUUmXXmpub7wBNGkCjBxpdUSWYJJERERE7k2dCpw4AVx5JXDWWVZH43dsbiMiIqKyxo0D9u41UwFMmmR1NJZgkkRERERl9elj1nKz24Hp04H9+62OyO+YJBEREVFZS5eemWW7qKhO1iYxSSIiIqKSsrJM7ZHNZrZttjpZm8QkiYiIiEqaNMk0s7mqg7VJTJKIiIiopBUrztQiOdlswPLl1sRjEU4BQERERCWtW2d1BAGBNUlEREREbjBJIiIiInKDSRIRERGRG0ySiIiIiNxgkkRERETkhqiq908qcghAutdPXFIcgMM+vkagqstlB+p2+ety2YG6XX6Wve6qy+X3V9kTVbVZ6Z0+SZL8QURWq2qK1XFYoS6XHajb5a/LZQfqdvlZ9rpZdqBul9/qsrO5jYiIiMgNJklEREREbgRzkjTN6gAsVJfLDtTt8tflsgN1u/wse91Vl8tvadmDtk8SERERkS8Fc00SERERkc8wSSIiIiJyIyCTJBG5SkS2icgOEXnCzeORIjLX8fjPIpLk8tiTjv3bRORKf8btDR6U/TER2SwiG0XkfyKS6PJYkYisd9wW+TfymvOg7KNE5JBLGe9xeWykiPzuuI30b+Te4UH533Ap+3YROebyWLC/9x+KyEER+a2cx0VE3nK8NhtF5HyXx4L6vfeg7MMdZf5VRJaLSHeXx9Ic+9eLyGr/Re0dHpT9MhHJcfndfsblsQo/L8HAg/KPcyn7b47PeRPHY8H+3rcRkW8d32ebROQRN8dY/7lX1YC6AQgFsBNAOwARADYA6FLqmLEA3nHcHwZgruN+F8fxkQDaOs4TanWZvFz2/gCiHffvd5bdsX3S6jL4uOyjAEx189wmAHY5fjZ23G9sdZm8Xf5Sxz8E4MPa8N474r8UwPkAfivn8cEAlgAQABcC+LkWvfeVlb2Ps0wArnaW3bGdBiDO6jL4sOyXAfjczf4qfV4C9VZZ+Usdey2Ab2rRe98CwPmO+w0AbHfzN9/yz30g1iT1BrBDVXepqg3AJwCuL3XM9QBmOO7PBzBARMSx/xNVzVfV3QB2OM4XLCotu6p+q6p5js2VAFr7OUZf8eR9L8+VAJaq6hFVPQpgKYCrfBSnr1S1/LcDmOOXyPxAVX8AcKSCQ64HMFONlQAaiUgL1IL3vrKyq+pyR9mA2vWZ9+R9L09N/l4EjCqWv7Z95rNUda3j/gkAWwC0KnWY5Z/7QEySWgHY47K9F2VfuOJjVLUQQA6Aph4+N5BVNf4/wGTZTlEislpEVorIDb4I0Ic8LfvNjmrX+SLSporPDWQel8HRxNoWwDcuu4P5vfdEea9PbXjvq6L0Z14BfCUia0RkjEUx+dpFIrJBRJaISFfHvjr1votINEwS8C+X3bXmvRfTZaYHgJ9LPWT55z7MFycl3xOREQBSAPRz2Z2oqpki0g7ANyLyq6rutCZCn/gPgDmqmi8i/wdTm3i5xTFZYRiA+apa5LKvtr/3dZ6I9IdJki5x2X2J431vDmCpiGx11E7UFmthfrdPishgAAsBnGNxTFa4FsBPqupa61Qr3nsRqQ+T/D2qqsetjqe0QKxJygTQxmW7tWOf22NEJAxAQwDZHj43kHkUv4gMBPA0gOtUNd+5X1UzHT93AfgOJjMPFpWWXVWzXcr7PoCenj43CFSlDMNQqto9yN97T5T3+tSG975SInIezO/89aqa7dzv8r4fBLAAwdW9oFKqelxVTzruLwYQLiJxqCPvu4uKPvNB+96LSDhMgpSqqv92c4j1n3srOmxVdIOp3doF05zg7JDXtdQxD6Bkx+1PHfe7omTH7V0Iro7bnpS9B0yHxXNK7W8MINJxPw7A7wiijowelr2Fy/0bAax03G8CYLfjNWjsuN/E6jJ5u/yO4zrBdNiU2vLeu5QjCeV34B2Ckh04f6kt770HZU+A6V/Zp9T+GAANXO4vB3CV1WXxctnPcv6uwyQBGY7fAY8+L8Fwq6j8jscbwvRbiqlN773jfZwJ4M0KjrH8cx9wzW2qWigiDwL4EmYEw4equklEXgCwWlUXAfgAwCwR2QHzyzPM8dxNIvIpgM0ACgE8oCWbJAKah2V/FUB9APNMX3VkqOp1ADoDeFdE7DA1hC+r6mZLClINHpb9YRG5Dua9PQIz2g2qekREJgFY5TjdC1qyWjrgeVh+wPyuf6KOvxQOQf3eA4CIzIEZyRQnInsBPAsgHABU9R0Ai2FGuuwAkAdgtOOxoH/vPSj7MzB9Lv/h+MwXqlkVPR7AAse+MAAfq+oXfi9ADXhQ9lsA3C8ihQBOARjm+N13+3mxoAg14kH5AfMP4Veqmuvy1KB/7wFcDOBOAL+KyHrHvqdg/ikImM89lyUhIvr/9u5dtYooCgPwv+wFGxvbYCEcvCA2PoKNIliJIoiNYBfBFxAEbbXwGey0txJrUUEQ8QkENSCkybI4A0bZVQyZQ873VXuurGmGnz3D2gADq/hPEgDA7IQkAIABIQkAYEBIAgAYEJIAAAaEJGBWVXWsqu5O4xNV9WLumgASLQCAmU3rNr3q7sXMpQD8ZeWaSQJr51GSjamh3Ockp7p7UVW3klzJsqPwySRPsuyufCPJdpJLU1O5jSRPkxzPsuHcne7+dPCPARw2PrcBc3uQ5Et3n01y/59jiyRXk1xI8jDJr+4+l+RtkpvTOc+T3Ovu80k2kzw7kKqBQ89MErDKXnf3VpKtqvqR5OW0/32S09MK4hfzZ5meZLl2I8B/E5KAVba9a7yza3sny/fXkSTfp1kogH3lcxswt60kR/dyYXf/TPK1qq4lSS2d2c/igPUlJAGz6u5vSd5U1Yckj/dwi+tJblfVuyQfk1zez/qA9aUFAADAgJkkAIABIQkAYEBIAgAYEJIAAAaEJACAASEJAGBASAIAGPgNjDWV4jmzoKwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "dztvJY33QeAM" }, "source": [ "The table below illustrates the absolute error and the signed error of the numerical method." ] }, { "cell_type": "code", "metadata": { "id": "Mb3KtEjTQeAM", "colab": { "base_uri": "https://localhost:8080/", "height": 357 }, "outputId": "28a671f5-41db-4806-fe04-2885b1e1177c" }, "source": [ "n=10\n", "d = {'time': t[0:n], 'Abysmal Kramer-Butler w': w[0:n],'Exact Error abs':np.abs(y[0:n]-w[0:n]),\n", " 'Exact Error':(y[0:n]-w[0:n])}\n", "df = pd.DataFrame(data=d)\n", "df" ], "execution_count": 9, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
timeAbysmal Kramer-Butler wExact Error absExact Error
00.0001.0000000.0000000.000000
10.1250.8899940.0000000.000000
20.2500.8675030.059902-0.059902
30.3750.7724910.022912-0.022912
40.5000.8231340.110072-0.110072
50.6250.7102970.014774-0.014774
60.7500.8612690.166535-0.166535
70.8750.6759860.0327380.032738
81.0000.9885920.252834-0.252834
91.1250.6319370.1423680.142368
\n", "
" ], "text/plain": [ " time Abysmal Kramer-Butler w Exact Error abs Exact Error\n", "0 0.000 1.000000 0.000000 0.000000\n", "1 0.125 0.889994 0.000000 0.000000\n", "2 0.250 0.867503 0.059902 -0.059902\n", "3 0.375 0.772491 0.022912 -0.022912\n", "4 0.500 0.823134 0.110072 -0.110072\n", "5 0.625 0.710297 0.014774 -0.014774\n", "6 0.750 0.861269 0.166535 -0.166535\n", "7 0.875 0.675986 0.032738 0.032738\n", "8 1.000 0.988592 0.252834 -0.252834\n", "9 1.125 0.631937 0.142368 0.142368" ] }, "metadata": {}, "execution_count": 9 } ] }, { "cell_type": "markdown", "metadata": { "id": "kDv2ORvuQeAN" }, "source": [ "# Theory\n", "## Consistent \n", "The Abysmal Kramer-Butler method does satisfy the consistency condition\n", "\\begin{equation}\\tau_{i}(h)=\\frac{y_{i+1}-y_{i-1}}{2h}-\\frac{1}{2}[4(f(t_i,y_i)-2f(t_{i-1},y_{i-1})] \\end{equation}\n", "As $h \\rightarrow 0$ \n", "\\begin{equation}\\frac{1}{2}[4(f(t_i,y_i)-2f(t_{i-1},y_{i-1})] \\rightarrow f(t_i,y_i).\\end{equation}\n", "While as $h \\rightarrow 0$ \n", "\\begin{equation}\\frac{y_{i+1}-y_{i-1}}{2h} \\rightarrow \\frac{y^{'}}{1}=\\frac{f(t_i,y_i)}{1}.\\end{equation}\n", "Hence as $h \\rightarrow 0$ $$\\frac{y_{i+1}-y_{i}}{h}-\\frac{1}{2}[4(f(t_i,y_i)-2f(t_{i-1},y_{i-1})]\\rightarrow f(t_i,y_i)-f(t_i,y_i)=0,\\end{equation}\n", "which means it is consistent.\n", "\n", "## Convergent \n", "The Abysmal Kramer-Butler method does satisfy the Lipschitz condition:\n", "\\begin{equation}F(t,w:h)-F(t,\\hat{w}:h)=\\frac{4}{2}f(t,w_i)-\\frac{2}{2}f(t-h,w_{i-1}))-(\\frac{4}{2}f(t,\\hat{w}_{i})-\\frac{2}{2}f(t-h,\\hat{w}_{i-1}))),\\end{equation}\n", "\\begin{equation}F(t,w:h)-F(t,\\hat{w}:h)=\\frac{4}{2}(f(t,w_i)-f(t,\\hat{w}_i))-\\frac{2}{2}(f(t-h,w_{i-1}))-f(t-h,\\hat{w}_{i-1}))),\\end{equation}\n", "\\begin{equation}F(t,w:h)-F(t,\\hat{w}:h)\\leq\\frac{4}{2}L|w_i-\\hat{w_i}|+\\frac{2}{2}L|w-\\hat{w}|\\leq \\frac{6}{2} L|w_i-\\hat{w_i}|,\\end{equation}\n", "This means it is internally convergent,\n", "\\begin{equation}|w_i-\\hat{w_i}|\\rightarrow 0\\end{equation}\n", "as $h \\rightarrow 0$.\n", "## Stability\n", "The Abysmal Kramer-Butler method does __not__ satisfy the stability condition.\n", "The characteristic equation of the \n", "\\begin{equation}w_{i+1} = w_{i-1} + h(4f(t_i,w_i)-2f(t_{i-1},w_{i-1})) \\end{equation}\n", "is\n", "\\begin{equation}\\lambda^2 = 1, \\end{equation}\n", "This has two roots $\\lambda=1$ and $\\lambda=-1$, hence the method is weakly stable." ] }, { "cell_type": "code", "metadata": { "id": "WeAzaYK2QeAO" }, "source": [ "" ], "execution_count": 9, "outputs": [] } ] }